A Broadscale Test of Host-Symbiont Cophylogeny Reveals Widespread Phylogenetic Congruence

Electronic Supplementary Material

Setups

Loading packages and custom functions

# loading packages
# devtools::install_github("thomasp85/patchwork")
pacman::p_load(tidyverse, # tidy family and related pacakges below
               kableExtra, 
               gridExtra, # may not use this
               purrr,
               magrittr, # extending piping
               pander,   # nice tables
               metafor,  # package for meta-analysis
               MCMCglmm,  # Bayeisan mixed model package
               ggbeeswarm, # making bee-swarm plots possible
               plotly,     # interactive plots using ggplot2
               MuMIn,  # multi-model inference
               lme4,   # lmm & glmm (models)
               broom.mixed, # getting estimates from lmer + glmer objects
               performance, # getting R2 from lmer + glmer objects
               png,         # reading png files
               grid,        # graphic layout manipulation
               patchwork    # putting ggplots together
               #lmerTest   # more functions for lme4
               #mi,      # missing data analysis
               #betareg   # dependance of the above
)

# getting functions 
source("../R/ESM_functions.R", chdir = TRUE)

Custom functions

We have 5 custom functions named : p_to_Zr(),I2(), R2(), get_est(), get_pred(), and cont_gen(), all of which are used later (see below for their functionality).

The Cophylogeny Dataset

Table of the dataset

Below is the dataset used for our meta-analysis followed by explanations of 24 variables original collected (not all variables were used for our analyses; variables which were neither ‘directly’ nor ‘indirectly’ used in our analyses are indicated by *).

Extended Data Table 1: the meta-analytic dataset of this study.

# getting the data and formating some variables (turning chraracter vectors
# to factors) read the difference between factr() and as.factor()
# <https://stackoverflow.com/questions/39279238/why-use-as-factor-instead-of-just-factor>
# full_data <- read.csv('../data/2019-04-04-source-data-dat.csv', na = 'NA',
# fileEncoding='UTF-8') %>% mutate_if(is.character, as.factor)
full_data <- read_csv("../data/2019-23-07-source-data-dat.csv", na = "NA") %>% 
    mutate_if(is.character, as.factor)

# making a scrollable table
kable(full_data, "html") %>% kable_styling("striped", position = "left") %>% 
    scroll_box(width = "100%", height = "500px")

authors	year	host_tax_broad	host_tax_fine	symbiont_tax_broad	symbiont_tax_fine	symbiont_euk	symbiosis	endo_or_ecto	mode_of_transmission_broad	mode_of_transmission_fine	symbiont	Visiting_symbiont?	host_tips_linked	host_tips_linked_corrected	host_genera	total_host_symbioint_links	host_range_link_ratio	host_range_taxonomic_breadth	symbiont_tips_linked	symbiont_genera	no_randomizations	p_value	method
Xu_et_al_2017	2017	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Buchnera	resident	20	20	8	20	1.00	1.00	20	1	1000	0.00100	TreeMap
Riess_et_al_2018	2018	Plant	Fungus	Microbe	Fungus	y	Parasite	Endo	horizontal	autonomous	Root fungus	resident	11	11	5	11	1.00	1.00	11	2	1000	0.33500	TreeMap
Badets_et_al_2011	2011	Vert	Tetrapod	Invert	Invert	y	Parasite	Endo	horizontal	autonomous	Monogenea	resident	17	17	13	17	1.00	1.00	17	4	1000	0.01000	TreeMap
Banks_et_al_2006	2006	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	Chewing lice	resident	18	18	6	30	2.00	1.60	15	2	1000	0.01000	TreeMap
Bochkov_et_al_2011	2011	Vert	Tetrapod	Invert	Invert	y	Parasite	Ecto	both	contact	Mites	resident	6	6	NA	9	1.00	1.00	9	6	100	0.01000	TreeMap
Charleston_&_Perkins_2003	2003	Vert	Tetrapod	Microbe	Protist	y	Parasite	Endo	horizontal	vector	Malaria	resident	9	9	1	9	1.00	1.00	9	1	100	0.03000	TreeMap
Charleston_&_Robertson_2002	2002	Vert	Tetrapod	Microbe	Virus	n	Parasite	Endo	horizontal	bodily fluid	Retrovirus	resident	12	12	12	12	1.00	1.00	12	1	1000	0.01500	TreeMap
Chauvatcharin_et_al_2006	2006	Microbe	Bacterium	Microbe	Virus	n	Parasite	Endo	both	NA	Bacteriophage	resident	14	14	7	15	1.07	1.00	14	1	1000	0.33100	TreeMap
Clark_et_al_2000	2000	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Buchnera	resident	17	17	4	17	1.00	1.00	17	1	1000	0.00100	TreeMap
Clayton_&_Johnson_2003	2003	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	Dove body lice	resident	13	13	6	14	1.08	1.15	13	1	10000	0.00060	TreeMap
Clayton_&_Johnson_2003	2003	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	Dove wing lice	resident	13	13	6	16	1.60	1.50	10	1	10000	0.15300	TreeMap
Cui_et_al_2012	2012	Vert	Tetrapod	Microbe	Virus	n	Parasite	Endo	horizontal	bodily fluid	Retrovirus	resident	7	7	3	7	1.00	1.00	7	7	10000	0.36600	TreeMap
Dabert_2001__	2001	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	Feather mite	resident	21	21	12	22	1.00	1.00	22	9	10000	0.00100	TreeMap
Deng_et_al_2013	2013	Invert	Invert	Invert	Invert	y	Parasite	Endo	horizontal	autonomous	Parasitic wasps	resident	7	7	4	10	1.00	1.43	10	2	999	0.53000	TreeMap
Desai_et_al_2010	2010	Microbe	Protist	Microbe	Bacterium	n	Mutualist	Ecto	both	contact	Devescovinid flagellates	resident	9	9	2	8	0.89	1.00	9	1	1000	0.00100	TreeMap
Desdevises_et_al_2002	2002	Vert	Fish	Invert	Invert	y	Parasite	Ecto	horizontal	autonomous	Monogenea	resident	14	14	11	39	1.95	1.45	20	2	999	0.31700	TreeMap
Downie_&_Gullan_2005	2005	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Primary endosymbiont Tremblaya	resident	21	21	16	21	1.00	1.00	21	1	1000	0.00100	TreeMap
Erpenbeck_et_al_2002	2002	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Sponge symbiont	resident	6	6	5	6	1.00	1.00	6	1	10000	0.02000	TreeMap
Etherington_et_al_2006	2006	Vert	Tetrapod	Microbe	Virus	n	Parasite	Endo	horizontal	contact	Calicivirus	resident	7	7	7	8	1.00	1.00	8	1	1000	0.01000	TreeMap
Farrell_1998	1998	Plant	Plant	Invert	Invert	y	Parasite	Ecto	horizontal	autonomous	Cerambicid beetle	resident	21	21	2	21	1.00	1.00	21	6	1000	0.07000	TreeMap
Gottschling_et_al_2011	2011	Vert	Tetrapod	Microbe	Virus	n	Parasite	Endo	horizontal	contact	Papillomavirus	resident	43	43	40	78	1.00	1.00	78	30	1000	0.00100	TreeMap
Hendricks_et_al_2013	2013	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	Quill mites	resident	19	19	17	20	1.25	1.25	16	1	1000	0.02100	TreeMap
Hosokawa_et_al_2006	2006	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Gamma proteobacteria	resident	7	7	3	7	1.00	1.00	7	1	1000	0.00100	TreeMap
Hugot_1999	1999	Vert	Tetrapod	Invert	Invert	y	Parasite	Endo	horizontal	trophic	Pinworms	resident	10	10	9	11	1.00	1.00	11	6	1000	0.00100	TreeMap
Hugot_et_al_2003	2003	Vert	Tetrapod	Microbe	Fungus	y	Parasite	Endo	horizontal	contact	Pneumocystis	resident	19	19	12	19	1.00	1.00	19	1	1000	0.00100	TreeMap
Huyse_&_Volckaert_2005	2005	Vert	Fish	Invert	Invert	y	Parasite	Ecto	horizontal	autonomous	Gyrodactylus flatworms	resident	8	8	3	22	1.29	1.24	17	1	100	0.02000	TreeMap
IkedaOhtsubo_&_Brune_2009	2009	Microbe	Protist	Microbe	Bacterium	n	Mutualist	Ecto	vertical	vertical	Endomicrobium	resident	11	11	1	11	1.00	1.00	11	1	1000	0.00100	TreeMap
Jackson_&_Charleston_1994	1994	Vert	Tetrapod	Microbe	Virus	n	Parasite	Endo	horizontal	contact	Lyssavirus	resident	10	10	8	10	1.00	1.00	10	1	100	0.19000	TreeMap
Jackson_&_Charleston_1994	1994	Vert	Tetrapod	Microbe	Virus	n	Parasite	Endo	horizontal	bodily fluid	Retrovirus	resident	10	10	7	10	1.00	1.00	10	1	100	0.01000	TreeMap
Jackson_&_Charleston_1994	1994	Vert	Tetrapod	Microbe	Virus	n	Parasite	Endo	horizontal	contact	Arenavirus	resident	11	11	7	12	1.00	1.00	12	1	100	0.05000	TreeMap
Jackson_&_Charleston_1994	1994	Vert	Bird	Microbe	Virus	n	Parasite	Endo	horizontal	bodily fluid	Retrovirus	resident	13	13	9	15	1.00	1.00	15	1	100	0.18000	TreeMap
Jackson_&_Charleston_1994	1994	Vert	Tetrapod	Microbe	Virus	n	Parasite	Endo	horizontal	contact	Hantavirus	resident	14	14	7	17	1.00	1.00	17	1	100	0.01000	TreeMap
Jeong_et_al_1999	1999	Plant	Plant	Microbe	Bacterium	n	Mutualist	Endo	horizontal	vector	Frankia	resident	12	12	12	19	1.06	1.00	18	1	1000	0.23000	TreeMap
Johnson_et_al_2002	2002	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	Bird louse	resident	25	25	23	25	1.32	1.42	19	19	1000	0.23000	TreeMap
Johnson_et_al_2003	2003	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	Dove wing lice	resident	28	28	15	31	1.48	1.33	21	1	100	0.03000	TreeMap
Johnson_et_al_2006	2006	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	Flamingo lice	resident	10	10	NA	11	1.00	1.00	11	11	1000	0.40200	TreeMap
Jousselin_et_al_2008	2009	Plant	Plant	Invert	Invert	y	Mutualist	Endo	horizontal	autonomous	Fig wasp	resident	15	15	1	15	1.07	1.07	14	6	10000	0.01000	TreeMap
Jousselin_et_al_2008	2009	Plant	Plant	Invert	Invert	y	Parasite	Endo	horizontal	autonomous	Fig wasp	resident	13	13	1	14	1.00	1.00	14	1	10000	0.01000	TreeMap
Jousselin_et_al_2008	2009	Plant	Plant	Invert	Invert	y	Parasite	Endo	horizontal	autonomous	Fig wasp	resident	13	13	1	13	1.00	1.00	13	1	10000	0.01000	TreeMap
Jousselin_et_al_2008	2009	Plant	Plant	Invert	Invert	y	Parasite	Endo	horizontal	autonomous	Fig wasp	resident	16	16	1	18	1.06	1.06	17	2	10000	0.01000	TreeMap
Jousselin_et_al_2009	2009	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Buchnera	resident	55	22	1	22	1.00	1.00	22	1	10000	0.00100	TreeMap
Kawaida_et_al_2013	2013	Invert	Invert	Plant	Plant	y	Mutualist	Endo	vertical	vertical	Green algae	resident	6	6	1	6	1.00	1.00	6	1	10000	0.00350	TreeMap
Kawakita_et_al_2004	2004	Plant	Plant	Invert	Invert	y	Mutualist	Ecto	horizontal	autonomous	Pollinating moth	visitor	18	18	1	18	1.00	1.00	18	1	999	0.01900	TreeMap
Kelley_&_Farrell_1998	1998	Plant	Plant	Invert	Invert	y	Parasite	Endo	horizontal	autonomous	Bark beetle	resident	41	41	1	89	6.85	1.92	13	1	100	0.28000	TreeMap
Kikuchi_et_al_2009	2009	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Gut symbiont	resident	14	14	5	14	1.00	1.00	14	1	1000	0.00100	TreeMap
Lanterbecq_et_al_2010	2010	Invert	Invert	Invert	Invert	y	Parasite	Endo/Ecto	horizontal	contact	Myzostomid worm	resident	16	16	12	16	1.00	1.00	16	5	5000	0.04000	TreeMap
Light_&_Hafner_2008	2008	Vert	Tetrapod	Invert	Invert	y	Parasite	Ecto	both	contact	Rodent sucking lice	resident	44	21	4	21	1.00	1.00	21	1	1000	0.00100	TreeMap
LimFong_et_al_2008	2008	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Endobugula bacteria	resident	5	5	1	5	1.00	1.00	5	2	1000	0.11000	TreeMap
Liu_et_al_2013	2013	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Buchnera	resident	37	37	1	37	1.00	1.00	37	1	1000	0.00100	TreeMap
Liu_et_al_2014	2014	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Buchnera	resident	27	27	3	29	1.00	1.00	29	1	1000	0.01000	TreeMap
LopezVaamonde_et_al_2001	2001	Invert	Invert	Invert	Invert	y	Parasite	Endo	horizontal	autonomous	Figwasp	resident	15	15	1	15	1.00	1.00	15	1	1000	0.00100	TreeMap
LopezVaamonde_et_al_2003	2003	Plant	Plant	Invert	Invert	y	Parasite	Endo	horizontal	autonomous	Leaf mining moth	resident	33	33	33	77	1.00	1.00	77	1	1000	0.21300	TreeMap
LopezVaamonde_et_al_2005	2005	Invert	Invert	Invert	Invert	y	Parasite	Endo	horizontal	autonomous	Parasitic wasp	resident	28	28	1	35	2.33	1.33	15	1	1000	0.24800	TreeMap
Martin_et_al_1999	1999	Vert	Tetrapod	Microbe	Virus	n	Parasite	Endo	horizontal	bodily fluid	Retrovirus	resident	38	38	38	48	1.00	1.00	48	1	1000	0.00100	TreeMap
Martin_et_al_2003	2003	Vert	Bird	Microbe	Virus	n	Parasite	Endo	horizontal	bodily fluid	Retrovirus	resident	14	14	14	16	1.00	1.00	16	1	100	0.01000	TreeMap
Martin_et_al_2003	2003	Vert	Tetrapod	Microbe	Virus	n	Parasite	Endo	horizontal	bodily fluid	Retrovirus	resident	9	9	9	9	1.00	1.00	9	1	100	0.05000	TreeMap
Martin_et_al_2003	2003	Vert	Tetrapod	Microbe	Virus	n	Parasite	Endo	horizontal	bodily fluid	Retrovirus	resident	17	17	17	23	1.00	1.00	23	1	100	0.21000	TreeMap
Mazzon_et_al_2010	2010	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Stammerula	resident	17	17	10	17	1.06	1.13	16	3	1000	0.00100	TreeMap
Morelli_&_Spicer_2007	2007	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	Bird nasal mite	resident	6	6	6	6	1.00	1.00	6	1	10000	0.00950	TreeMap
Muniz_et_al_2013	2013	Vert	Tetrapod	Microbe	Virus	n	Parasite	Endo	horizontal	bodily fluid	Retrovirus	resident	17	17	11	17	1.00	1.00	17	1	10000	0.00001	TreeMap
Musser_et_al_2010	2010	Vert	Mammal	Invert	Invert	y	Parasite	Ecto	both	contact	Squirrel sucking lice	resident	6	6	1	6	1.00	1.00	6	6	100	0.30000	TreeMap
Pagan_et_al_2010	2010	Plant	Plant	Microbe	Virus	n	Parasite	Endo	horizontal	contact	Tobamovirus	resident	10	10	10	13	1.00	1.00	13	1	1000	0.24000	TreeMap
Page_1996	1996	Vert	Tetrapod	Invert	Invert	y	Parasite	Ecto	both	contact	Gopher chewing lice	resident	15	15	6	17	1.00	1.00	17	2	1000	0.00100	TreeMap
Page_et_al_1998	1998	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	Bird lice	resident	7	7	1	8	1.00	1.00	8	1	100	0.01000	TreeMap
Page_et_al_2004	2004	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	Seabird lice	resident	12	12	4	12	1.00	1.00	12	1	1000	0.00100	TreeMap
Page_et_al_2004	2004	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	Seabird lice	resident	11	11	5	13	1.00	1.00	13	4	100	0.25000	TreeMap
Page_et_al_2004	2004	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	Seabird lice	resident	13	13	5	14	1.00	1.00	14	1	100	0.46000	TreeMap
Page_et_al_2004	2004	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	Seabird lice	resident	9	9	3	9	1.00	1.00	9	1	100	0.36000	TreeMap
Paterson_&_Poulin_1999	1999	Vert	Fish	Invert	Invert	y	Parasite	Ecto	horizontal	autonomous	Parasitic copepod	resident	8	8	8	12	1.20	2.20	10	1	100	0.01000	TreeMap
Paterson_et_al_2000	2000	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	Seabird lice	resident	11	11	5	14	1.00	1.00	14	5	100	0.01000	TreeMap
Percy_et_al_2004	2004	Plant	Plant	Invert	Invert	y	Parasite	Ecto	both	autonomous	Psyllid	resident	35	35	8	56	1.22	1.09	46	4	1000	0.00500	TreeMap
PerezLosada_et_al_2006	2006	Vert	Tetrapod	Microbe	Virus	n	Parasite	Endo	horizontal	NA	Polyomavirus	resident	9	9	9	11	1.00	1.00	11	1	100	0.01000	TreeMap
Perlman_et_al_2003	2003	Invert	Invert	Invert	Invert	y	Parasite	Endo	horizontal	autonomous	Nematode	resident	16	16	1	17	1.89	1.33	9	1	1000	0.14000	TreeMap
Quek_et_al_2004	2004	Plant	Plant	Invert	Invert	y	Mutualist	Ecto	horizontal	autonomous	Ants	resident	11	11	1	23	2.30	1.60	10	1	1000	0.99900	TreeMap
Ramsden_et_al_2008	2008	Vert	Tetrapod	Microbe	Virus	n	Parasite	Endo	horizontal	contact	Hantavirus	resident	33	33	20	38	1.00	1.00	38	1	1000	0.99900	TreeMap
Reed_et_al_2007	2007	Vert	Tetrapod	Invert	Invert	y	Parasite	Ecto	horizontal	contact	Anthropoid lice	resident	4	4	4	6	1.20	1.00	5	4	1000	0.05000	TreeMap
Refregier_et_al_2008	2008	Plant	Plant	Microbe	Fungus	y	Parasite	Endo	horizontal	vector	Anther smut fungi	resident	18	18	7	20	1.25	1.25	16	1	3000	0.50000	TreeMap
Shoemaker_et_al_2002	2002	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	both	NA	Wolbachia	resident	20	20	7	23	1.00	1.00	23	1	100	0.47000	TreeMap
Simkova_et_al_2013	2013	Vert	Fish	Invert	Invert	y	Parasite	Ecto	horizontal	autonomous	Monogenea	resident	5	5	1	21	1.00	1.00	21	1	999	0.04500	TreeMap
Six_&_Paine_1999	1999	Microbe	Fungus	Invert	Invert	y	Mutualist	Ecto	horizontal	vector	Mycangial fungi	resident	6	6	1	6	1.00	1.00	6	3	1000	0.03100	TreeMap
Skerikova_et_al_2001	2001	Vert	Fish	Invert	Invert	y	Parasite	Endo	horizontal	trophic	Cestode	resident	7	7	7	7	1.00	1.00	7	1	10000	0.40000	TreeMap
Smith_et_al_2008b	2008	Vert	Tetrapod	Invert	Invert	y	Parasite	Ecto	both	contact	Rodent lice	resident	20	20	3	20	1.00	1.00	20	14	100	0.05000	TreeMap
Sorenson_et_al_2004	2004	Vert	Bird	Vert	Bird	y	Parasite	Ecto	horizontal	autonomous	Brood parasitic finch	resident	33	33	10	34	1.62	1.43	21	1	100	0.15000	TreeMap
Subbotin_et_al_2004	2004	Plant	Plant	Invert	Invert	y	Parasite	Endo	horizontal	autonomous	Nematode	resident	16	16	16	21	1.00	1.00	21	4	1000	0.00100	TreeMap
Switzer_et_al_2005	2005	Vert	Tetrapod	Microbe	Virus	n	Parasite	Endo	horizontal	bodily fluid	Retrovirus	resident	46	46	17	51	1.00	1.00	51	1	10000	0.00700	TreeMap
Urban_&_Cryan_2012	2012	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	autonomous	Fulgorid planthopper microbe	resident	40	40	38	40	1.00	1.00	40	1	1000	0.00100	TreeMap
Urban_&_Cryan_2012	2012	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	autonomous	Fulgorid planthopper microbe	resident	30	30	29	30	1.00	1.00	30	1	1000	0.00100	TreeMap
Vanhove_et_al_2015	2015	Vert	Fish	Invert	Invert	y	Parasite	Endo	horizontal	autonomous	Monogenea	resident	19	19	10	28	1.00	1.00	28	1	10000	0.04210	TreeMap
Weckstein_2004	2004	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	Toucan chewing lice	resident	11	11	1	11	2.20	1.80	5	1	10000	0.89000	TreeMap
Weiblen_&_Bush_2002	2002	Plant	Plant	Invert	Invert	y	Mutualist	Endo	horizontal	autonomous	Fig wasp	resident	19	19	1	19	1.00	1.00	19	1	10000	0.01950	TreeMap
Weiblen_&_Bush_2002	2002	Plant	Plant	Invert	Invert	y	Parasite	Endo	horizontal	autonomous	Fig wasp	resident	12	12	1	18	1.00	1.00	18	1	10000	0.12150	TreeMap
Wu_et_al_2008	2008	Plant	Plant	Microbe	Virus	n	Parasite	Endo	horizontal	vector	Mastrevirus	resident	8	8	8	10	1.00	1.00	10	1	100	0.01000	TreeMap
Yan_et_al_2011	2011	Vert	Fish	Microbe	Virus	n	Parasite	Endo	horizontal	contact	Lymphocystis	resident	8	8	8	15	1.00	1.00	15	1	999	0.98900	TreeMap
Riess_et_al_2018	2018	Plant	Fungus	Microbe	Fungus	y	Parasite	Endo	horizontal	autonomous	Root fungus	resident	11	11	5	11	1.00	1.00	11	2	9999	0.26400	Parafit
Souza_et_al_2018	2018	Vert	Mammal	Microbe	Virus	n	Parasite	Endo	both	bodily fluid	Primate hepadnaviruses	resident	8	8	8	19	1.06	1.00	18	1	999	0.00500	Parafit
Ramasindrazana_et_al_2017	2017	Vert	Mammal	Invert	Invert	y	Parasite	Ecto	both	contact	Bat flies (Nycteribiidae)	resident	15	15	6	26	2.89	1.78	9	5	999	0.00100	Parafit
Li_et_al_2017	2017	Plant	Plant	Microbe	Fungus	y	Parasite	Endo	horizontal	vector	Grass fungus (Tranzscheliella)	Resident	12	12	9	12	1.71	2.00	7	1	999	0.50505	Parafit
Sweet_et_al_2018b	2018	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	body lice (Physconelloides)	resident	11	11	4	13	1.86	1.43	7	1	100000	0.00500	Parafit
Sweet_&_Johnson_2018	2018	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	wing lice (Columbicola)	resident	13	13	4	14	2.80	1.80	5	1	100000	0.00500	Parafit
Arab_et_al_2019	2019	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Blattabacterium of cockroaches	resident	55	55	52	55	1.00	1.00	55	1	999	0.00100	Parafit
Sweet_et_al_2018a	2018	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	Brueelia bird lice	resident	259	259	163	283	1.63	1.61	174	11	9999	0.00010	Parafit
Hewitt_et_al_2019	2019	Vert	Fish	Invert	Invert	y	Parasite	Ecto	horizontal	autonomous	Unionid mussels	resident	178	178	79	495	7.17	3.30	69	35	999	0.00100	Parafit
Latinne_et_al_2018	2018	Vert	Mammal	Microbe	Funus	y	Parasite	Endo	horizontal	enviromental	Pneumocystis of rodents	resident	15	15	6	19	3.17	2.00	6	1	999	0.00900	Parafit
Graca_et_al_2018	2018	Vert	Fish	Invert	Invert	y	Parasite	Endo	horizontal	autonomous	Monogenea	resident	9	9	7	18	1.29	1.29	14	1	10000	0.00010	Parafit
Sweet_et_al_2017	2017	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	Phabine bird lice	resident	12	12	5	15	1.15	1.23	13	3	999	0.06900	Parafit
Zhang_et_al_2017	2017	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Nardonella of Weevils	resident	44	44	24	44	1.00	1.00	44	1	100000	0.00001	Parafit
Megia-palma_et_al_2018	2018	Vert	Reptile	Microbe	Protist	y	Parasite	Endo	horizontal	vector	Lizard Schellackia Apicomplexa	resident	16	16	8	23	1.53	1.27	15	8	999	0.00400	Parafit
Dona_2018	2018	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	Trouessartia Feather mite	resident	14	14	13	15	1.00	1.00	15	1	100000	0.01000	Parafit
Dona_2018	2018	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	Proctophyllodes Feather mite	resident	42	42	29	44	1.00	1.00	44	1	100000	0.01000	Parafit
Xu_et_el_2017	2017	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Buchner	resident	20	20	8	20	1.00	1.00	20	1	999	0.00100	Parafit
Chen_et_al_2017	2017	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Buchnera	resident	50	50	11	50	1.00	1.00	50	1	999	0.00100	Parafit
Patra_et_al_2018	2018	Vert	Fish	Invert	Invert	y	Parasite	Endo	horizontal	environmental	Myxozoa	resident	24	24	24	31	1.00	1.00	31	1	999	0.00100	Parafit
Matthews_et_al_2018	2018	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	Amerodectes feather mites	resident	12	12	6	33	1.00	0.58	33	1	99900	0.00400	Parafit
Catanach_et_al_2018	2018	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	Colpocephalum	resident	54	28	5	44	1.07	0.93	41	30	999	0.00100	Parafit
Ballinger_et_al_2018	2018	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Spiroplasma	resident	11	11	1	12	1.00	1.00	12	1	999	0.07000	Parafit
Li_et_al_2018	2018	Vert	Vert	Invert	Invert	y	Parasite	Endo	horizontal	trophic	Ascarid worms	resident	68	68	34	129	1.45	NA	89	89	100000	0.00100	Parafit
Holzer_et_al_2018	2018	Invert	Invert	Invert	Invert	y	Parasite	Endo	horizontal	environmental	Myxozoa	resident	23	23	22	39	1.00	1.00	39	21	1000	0.00100	Parafit
Holzer_et_al_2018	2018	Vert	Fish	Invert	Invert	y	Parasite	Endo	horizontal	environmental	Myxozoa FW (OIM)	resident	69	69	62	101	1.00	1.00	101	15	1000	0.00100	Parafit
Holzer_et_al_2018	2018	Vert	Fish	Invert	Invert	y	Parasite	Endo	horizontal	environmental	Myxozoa (PIM)	resident	69	69	58	75	1.00	1.00	75	21	1000	0.00100	Parafit
Carneiro_et_al_2018	2018	Vert	Mammal	Microbe	Virus	n	Parasite	Endo	horizontal	bodily fluid	Hepatovirus	resident	26	26	21	26	1.00	1.00	26	1	100000	0.01000	Parafit
Jesovnik_et_al_2017	2017	Invert	Invert	Microbe	Fungus	y	Mutualist	Ecto	both	environmental	Ant fungus	resident	32	6	1	6	1.00	1.00	6	1	999	0.02900	Parafit
Singh_et_al_2017	2017	Microbe	Fungus	Plant	Plant	y	Mutualist	Endo	vertical	vertical	Lichen algae	resident	23	17	1	25	1.25	1.15	20	1	9999	0.00020	Parafit
Endara_et_al_2018	2018	Plant	Plant	Invert	Invert	y	Parasite	Endo	horizontal	autonomous	Sawfly	resident	44	30	1	45	1.18	1.03	38	NA	9999	0.01500	Parafit
Endara_et_al_2017	2017	Plant	Plant	Invert	Invert	y	Parasite	Ecto	horizontal	autonomous	Gelechioidea moths	visitor	18	18	1	29	1.00	2.17	29	NA	100	0.70000	Parafit
Endara_et_al_2017	2017	Plant	Plant	Invert	Invert	y	Parasite	Ecto	horizontal	autonomous	Riodinidae moths	visitor	10	10	1	18	1.50	1.25	12	NA	100	0.90000	Parafit
Endara_et_al_2017	2017	Plant	Plant	Invert	Invert	y	Parasite	Ecto	horizontal	autonomous	Erebidae moths	visitor	14	14	1	18	1.20	1.13	15	NA	100	0.74000	Parafit
Lauber_et_al_2017	2017	Vert	Vert	Microbe	Virus	n	Parasite	Endo	horizontal	autonomous	Nackednaviruses and Hepadnaviruses	resident	30	30	28	30	0.88	1.00	34	8	10000	0.00010	Parafit
Althoff_et_al_2012	2012	Plant	Plant	Invert	Invert	y	Mutualist	Ecto	horizontal	autonomous	Yucca moth	visitor	24	24	1	40	2.00	1.35	20	1	1000	0.00100	Parafit
Althoff_et_al_2012	2012	Plant	Plant	Invert	Invert	y	Mutualist	Ecto	horizontal	autonomous	Yucca moth	visitor	24	24	1	38	2.24	1.53	17	1	1000	0.00100	Parafit
Banks_et_al_2006	2006	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	Penguin chewing lice	resident	18	18	6	30	2.00	1.60	15	2	10000	0.00100	Parafit
Bayerlova_2009	2009	Vert	Mammal	Microbe	Virus	n	Parasite	Endo	horizontal	contact	Arenavirus	resident	21	21	14	31	1.55	1.45	20	1	9999	0.00040	Parafit
Bellec_et_al_2014	2014	Plant	Plant	Microbe	Virus	n	Parasite	Endo	horizontal	contact	Prasinovirus	resident	22	22	3	133	2.61	1.65	51	1	999	0.00100	Parafit
Bruyndonckxx_et_al_2009	2009	Vert	Mammal	Invert	Invert	y	Parasite	Ecto	both	contact	Bat mites	resident	20	20	7	21	1.91	2.27	11	2	9999	0.00300	Parafit
Caraguel_et_al__2007	2007	Microbe	Amoeba	Microbe	Protist	y	Mutualist	Endo	vertical	vertical	Prokinetoplastid endosymbiont	resident	6	6	1	6	1.00	1.00	6	1	9999	0.00100	Parafit
Choi_&_Thines_2015	2015	Plant	Plant	Microbe	Fungus	y	Parasite	Endo	horizontal	autonomous	Downy mildew	resident	63	63	28	63	1.00	1.00	63	3	999	0.00100	Parafit
Conord_et_al_2008	2008	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Dryophthoridae endosymbionts	resident	14	14	10	14	1.00	1.00	14	1	999	0.00900	Parafit
Cornuault_et_al_2012	2012	Vert	Bird	Microbe	Protist	y	Parasite	Endo	horizontal	vector	Leucocytozoon	resident	8	8	1	23	1.28	1.17	18	1	9999	0.03500	Parafit
Cruaud_et_al_2012	2012	Plant	Plant	Invert	Invert	y	Mutualist	Endo	horizontal	autonomous	Fig wasp	resident	200	200	1	200	1.00	1.00	200	20	9999	0.01000	Parafit
Cui_et_al_2014	2014	Vert	Bird	Microbe	Virus	n	Parasite	Endo	horizontal	bodily fluid	Hepadnavirus	resident	46	46	46	NA	NA	NA	61	NA	99999	0.23300	Parafit
Deng_et_al_2013	2013	Invert	Invert	Invert	Invert	y	Parasite	Endo	horizontal	autonomous	Parasitic wasp	resident	7	7	4	10	1.00	1.00	10	2	999	0.01602	Parafit
Desdevises_et_al_2002	2002	Vert	Fish	Invert	Invert	y	Parasite	Ecto	horizontal	autonomous	Monogenea	resident	14	14	11	39	2.00	1.65	20	2	999	0.26000	Parafit
Dhami_et_al_2013	2013	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Hoataupuhia symbiont	resident	42	42	NA	42	1.00	1.00	42	4	999	0.00100	Parafit
Dowie_et_al_2016	2016	Microbe	Fungus	Plant	Plant	y	Parasite	Endo	horizontal	environmental	Parasitic plant	resident	4	4	1	9	2.25	1.50	4	1	9999	0.00100	Parafit
Du_Toit_et_al_2013	2013	Vert	Mammal	Invert	Invert	y	Parasite	Ecto	both	contact	Rodent lice	resident	4	4	1	14	1.17	1.17	12	1	10000	0.88000	Parafit
FerrerParis_et_al_2013	2013	Plant	Plant	Invert	Invert	y	Parasite	Ecto	horizontal	autonomous	Butterflies	visitor	64	64	NA	112	2.73	5.20	41	NA	999	0.15700	Parafit
FraijaFernandez_et_al_2016	2016	Vert	Mammal	Invert	Invert	y	Parasite	Endo	horizontal	trophic	Digenea	resident	31	31	24	50	5.56	3.67	9	6	999	0.00100	Parafit
Garamszegi_2009	2009	Vert	Mammal	Microbe	Protist	y	Parasite	Endo	horizontal	vector	Primate malaria	resident	23	23	23	43	2.39	2.56	18	1	1000	0.00100	Parafit
Garcia_&_Hayman_2016	2016	Vert	Vert	Microbe	Protist	y	Parasite	Endo	horizontal	trophic	Cryptosporidium	resident	22	22	22	36	1.33	1.96	27	1	999	0.01000	Parafit
Gavotte_et_al_2007	2007	Microbe	Bacterium	Microbe	Virus	n	Parasite	Endo	both	NA	Bacteriophage	resident	33	33	1	51	0.93	1.00	55	1	10000	0.13190	Parafit
Goker_et_al_2011	2011	Microbe	Fungus	Microbe	Virus	n	Parasite	Endo	NA	NA	Mycovirus	resident	8	8	8	8	1.00	1.00	8	1	9999	0.09780	Parafit
Gomard_et_al_2016	2016	Vert	Mammal	Microbe	Bacterium	n	Parasite	Endo	horizontal	NA	Leptospira	resident	12	12	11	26	1.00	1.00	26	1	999	0.09000	Parafit
Gottschling_et_al_2011	2011	Vert	Mammal	Microbe	Virus	n	Parasite	Endo	horizontal	contact	Papillomavirus	resident	43	43	40	76	0.97	1.00	78	30	9999	0.00010	Parafit
Hall_et_al__2016	2016	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Carsonella primary endosymbiont	resident	37	37	18	37	1.00	1.00	37	1	10000	0.00100	Parafit
Hall_et_al__2016	2016	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Arsenophonus secondary endosymbiont	resident	20	20	9	20	1.00	1.00	20	1	10000	0.38700	Parafit
Hammer_et_al_2010	2010	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	Petrel chewing lice	resident	23	23	6	23	1.10	1.00	21	1	999	0.00100	Parafit
Hammerlinck_et_al_2016	2016	Invert	Invert	Invert	Invert	y	Parasite	Endo	horizontal	autonomous	Parasitic wasp	resident	11	11	2	11	1.38	1.13	8	2	999	0.12900	Parafit
Hammerlinck_et_al_2016	2016	Invert	Invert	Invert	Invert	y	Parasite	Endo	horizontal	autonomous	Parasitic wasp	resident	6	6	2	6	1.50	1.25	4	2	999	0.19500	Parafit
Hammerlinck_et_al_2016	2016	Invert	Invert	Invert	Invert	y	Parasite	Endo	horizontal	autonomous	Parasitic wasp	resident	7	7	2	7	1.00	1.00	7	2	999	0.04200	Parafit
Hembry_et_al_2013	2013	Plant	Plant	Invert	Invert	y	Mutualist	Ecto	horizontal	autonomous	Pollinating moths	visitor	37	37	1	35	1.00	1.00	35	1	999	0.01700	Parafit
Herrera__et_al_2016	2016	Microbe	Fungus	Microbe	Fungus	y	Parasite	Endo	horizontal	environmental	Mycoparasite Cosmospora	resident	13	13	9	13	1.00	1.00	13	3	999	0.00500	Parafit
Hoglund_et_al_2003	2003	Vert	Vert	Invert	Invert	y	Parasite	Endo	horizontal	trophic	Lungworm	resident	8	8	8	10	2.00	2.40	5	1	10000	0.09800	Parafit
Hughes_et_al_2007	2007	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	Pelican lice	resident	18	18	6	18	1.06	1.00	17	1	999	0.00010	Parafit
Huyse_&_Volckaert_2005	2005	Vert	Fish	Invert	Invert	y	Parasite	Ecto	horizontal	autonomous	Monogenea	resident	8	8	5	22	1.29	1.29	17	1	999	0.09500	Parafit
Irwin_et_al_2012	2012	Vert	Mammal	Microbe	Virus	n	Parasite	Endo	horizontal	contact	Arenavirus	resident	21	21	14	31	1.55	1.85	20	NA	9999	0.00040	Parafit
Jenkins_et_al_2012	2012	Vert	Bird	Microbe	Protist	y	Parasite	Endo	horizontal	autonomous	Leucocytozoon	resident	52	52	40	138	1.45	NA	95	1	999	0.00100	Parafit
Jousselin_et_al_2008	2008	Plant	Plant	Invert	Invert	y	Mutualist	Endo	horizontal	autonomous	Fig wasp	resident	15	15	1	15	1.07	1.07	14	6	9999	0.00300	Parafit
Jousselin_et_al_2008	2008	Plant	Plant	Invert	Invert	y	Parasite	Endo	horizontal	autonomous	Fig wasp	resident	13	13	1	13	1.00	1.00	13	1	9999	0.02000	Parafit
Jousselin_et_al_2008	2008	Plant	Plant	Invert	Invert	y	Parasite	Endo	horizontal	autonomous	Fig wasp	resident	16	16	1	18	1.06	1.06	17	2	9999	0.00100	Parafit
Jousselin_et_al_2008	2008	Plant	Plant	Invert	Invert	y	Parasite	Endo	horizontal	autonomous	Fig wasp	resident	13	13	1	14	1.00	1.00	14	1	9999	0.00300	Parafit
Jousselin_et_al_2009	2009	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Buchnera	resident	55	22	1	22	1.00	1.00	22	1	9999	0.00100	Parafit
Kaltenpoth_et_al_2014	2014	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Beewolf bacterial symbiont	resident	39	39	2	41	1.00	1.00	41	1	1000	0.00100	Parafit
Kawakita_&_Kato_2009	2009	Plant	Plant	Invert	Invert	y	Mutualist	Ecto	horizontal	autonomous	Yucca moth	visitor	10	10	7	10	1.00	1.00	10	1	100	0.39000	Parafit
Kawakita_et_al_2004	2004	Plant	Plant	Invert	Invert	y	Mutualist	Ecto	horizontal	autonomous	Pollinating moth	visitor	18	18	1	18	1.00	1.00	18	1	999	0.00500	Parafit
Kawazoe_et_al_2008	2008	Invert	Invert	Invert	Invert	y	Mutualist	Ecto	vertical	vertical	Bee mite	resident	4	4	1	5	1.25	1.25	4	1	9999	0.00010	Parafit
Kolsch_&__Pedersen_2010	2010	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Reed beetle bacterial endosymbiont	resident	41	41	7	41	1.00	1.00	41	1	9999	0.00100	Parafit
Krasnov_&_Shenbrot_2013	2013	Vert	Mammal	Invert	Invert	y	Parasite	Ecto	horizontal	contact	Jerboa fleas	resident	21	21	8	62	3.26	4.37	19	7	999	0.16000	Parafit
Krumbholz_et_al_2009	2009	Vert	Vert	Microbe	Virus	n	Parasite	Endo	horizontal	NA	Polyomavirus	resident	13	13	13	18	1.00	1.00	18	1	999999	0.49460	Parafit
Ku_&_Hu_2014	2014	Plant	Plant	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Burkholderia plant symbiont	resident	11	11	1	11	1.00	1.00	11	1	9999	0.02000	Parafit
Lanterbecq_et_al_2010	2010	Invert	Invert	Invert	Invert	y	Parasite	Endo/Ecto	horizontal	contact	Myzostomid worm	resident	16	16	12	16	1.00	1.00	16	5	1000	0.01300	Parafit
Lauron_et_al_2015	2015	Vert	Bird	Microbe	Protist	y	Parasite	Endo	horizontal	vector	Avian malaria	resident	18	18	8	83	1.80	NA	46	1	1000	0.66000	Parafit
Lei_&_Olival_2014	2014	Vert	Mammal	Microbe	Bacterium	n	Parasite	Endo	horizontal	environmental	New World bat Bartonella	resident	14	14	11	38	1.00	1.00	38	1	999	0.00100	Parafit
Lei_&_Olival_2014	2014	Vert	Mammal	Microbe	Bacterium	n	Parasite	Endo	horizontal	environmental	New World rodent Bartonella	resident	4	4	4	20	1.00	1.00	20	1	999	0.00100	Parafit
Lei_&_Olival_2014	2014	Vert	Mammal	Microbe	Bacterium	n	Parasite	Endo	horizontal	environmental	New World bat Leptospira	resident	14	14	12	19	1.00	1.00	19	1	999	0.85800	Parafit
Lei_&_Olival_2014	2014	Vert	Mammal	Microbe	Bacterium	n	Parasite	Endo	horizontal	environmental	Old World bat Bartonella	resident	9	9	8	13	1.00	1.00	13	1	999	0.02900	Parafit
Lei_&_Olival_2014	2014	Vert	Mammal	Microbe	Bacterium	n	Parasite	Endo	horizontal	environmental	Old World rodent Bartonella	resident	35	35	22	119	1.09	1.17	109	1	999	0.00010	Parafit
Lei_&_Olival_2014	2014	Vert	Mammal	Microbe	Bacterium	n	Parasite	Endo	horizontal	environmental	Old World bat Leptospira	resident	6	6	5	7	1.00	1.00	7	1	999	0.75870	Parafit
LewisRogers_&_Crandall_2009	2009	Vert	Tetrapod	Microbe	Virus	n	Parasite	Endo	horizontal	contact	Picornaviridae	resident	6	6	NA	27	1.00	1.00	27	11	999	0.47000	Parafit
Li_et_al_2017	2017	Plant	Plant	Microbe	Fungus	y	Parasite	Endo	horizontal	autonomous	Smut fungi	resident	12	12	9	12	1.71	1.57	7	1	999	0.50505	Parafit
Light_&_Hafner_2008	2008	Vert	Mammal	Invert	Invert	y	Parasite	Ecto	both	contact	Rodent sucking lice	resident	44	21	4	21	1.00	1.00	21	1	999	0.00100	Parafit
Liu_et_al_2016	2016	Plant	Plant	Microbe	Fungus	y	Parasite	Endo	NA	environmental	Tree foliar fungi	resident	13	13	10	44	1.57	2.71	28	28	9999	0.02510	Parafit
Liu_et_al_2016	2016	Plant	Plant	Microbe	Fungus	y	Parasite	Endo	NA	environmental	Subtropical tree soil fungi	resident	19	19	16	76	3.30	3.22	23	23	9999	0.02030	Parafit
Maneesakorn_et_al_2011	2011	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Photorhabdus symbiotic bacteria in nematodes	resident	12	12	1	12	1.00	1.00	12	1	999	0.00100	Parafit
Martinez_et_al_2011	2011	Invert	Invert	Microbe	Protist	y	Parasite	Endo	horizontal	vector	Avian malaria	resident	16	16	1	22	2.44	1.56	9	1	999	0.07500	Parafit
Martinez_San_udo_&_Girolami_2009	2009	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Tephritid symbiotic bacteria	resident	19	19	10	19	1.12	1.18	17	5	999	0.00300	Parafit
Mattiucci_&_nascetti_2008	2008	Vert	Mammal	Invert	Invert	y	Parasite	Endo	horizontal	trophic	Nematode	resident	7	7	15	12	1.33	1.78	9	1	100	0.05000	Parafit
Mazzon_et_al_2010	2010	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Stammerula bacterial symbiont	resident	17	17	10	17	1.06	1.00	16	3	999	0.00700	Parafit
McFrederick_&_Taylor_2013	2013	Invert	Invert	Invert	Invert	y	NA	Ecto	vertical	vertical	Nematode	resident	7	7	3	7	1.00	1.00	7	1	10000	0.00600	Parafit
McKee_et_al_2016	2016	Vert	Mammal	Microbe	Bacterium	n	Parasite	Endo	horizontal	vector	Bat bartonella bacteria	resident	66	66	42	184	1.06	1.09	173	1	10000	0.00010	Parafit
McLeish_&_Van_Noort_2012	2012	Plant	Plant	Invert	Invert	y	Mutualist	Endo	horizontal	autonomous	Fig wasp	resident	26	26	1	65	1.00	1.00	65	6	9999	0.18000	Parafit
Mehdiabadi_et_al_2012	2012	Invert	Invert	Microbe	Fungus	y	Mutualist	Ecto	both	environmental	Ant fungus	resident	99	11	1	11	1.00	1.00	11	1	999	0.00100	Parafit
Mendlova_et_al_2012	2012	Vert	Fish	Invert	Invert	y	Parasite	Ecto	horizontal	autonomous	Monogenea	resident	6	6	5	34	1.17	1.21	29	2	999	0.02400	Parafit
Merville_et_al._2013	2013	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Curculioniphilus buchneri primary endosymbiont of weevils	resident	9	9	1	9	1.00	1.00	9	1	999	0.00350	Parafit
Millanes_et_al_2014	2014	Microbe	Fungus	Microbe	Fungus	y	Parasite	Endo	both	environmental	Biatoropsis fungus	resident	16	16	2	16	1.00	1.00	16	1	999	0.23600	Parafit
Miyaki_et_al_2016	2016	Vert	Fish	Microbe	Bacterium	n	Mutualist	Endo	horizontal	environmental	Giant bacteria Epulopiscium	resident	8	8	4	54	3.18	2.24	17	1	999	0.00100	Parafit
Mondo_et_al_2012	2012	Microbe	Fungus	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Glomeribacter gigasporarum fungal symbiont	resident	55	5	2	5	1.00	1.00	5	1	10000	0.00010	Parafit
Nouioui_et_al_2014	2014	Plant	Plant	Microbe	Bacterium	n	Mutualist	Endo	horizontal	vector	Frankia	resident	9	9	1	20	1.00	1.00	20	1	9999	0.33000	Parafit
Pellissier_et_al_2013	2013	Plant	Plant	Invert	Invert	y	Parasite	Ecto	horizontal	autonomous	Butterflies	visitor	104	104	NA	NA	NA	NA	97	NA	10000	0.00010	Parafit
Perkins_2010	2010	Vert	Fish	Invert	Invert	y	Parasite	Endo	horizontal	autonomous	Monogenea	resident	61	61	50	75	1.00	1.00	75	NA	9999	0.92600	Parafit
Peterson_et_al_2010	2010	Plant	Plant	Microbe	Fungus	y	Parasite	Endo	horizontal	environmental	Beech fungus	resident	11	11	1	25	2.08	1.58	12	1	9999	0.00010	Parafit
Polme_et_al_2014.pdf	2014	Plant	Plant	Microbe	Bacterium	n	Mutualist	Endo	horizontal	vector	Frankia	resident	22	22	1	NA	NA	NA	43	1	999	0.00100	Parafit
Quek_et_al_2004	2004	Plant	Plant	Invert	Invert	y	Mutualist	Ecto	horizontal	autonomous	Ant	resident	10	10	1	23	2.30	1.70	10	1	999	0.77900	Parafit
Ricklefs_et_al_2004	2004	Vert	Bird	Microbe	Protist	y	Parasite	Endo	horizontal	vector	Avian malaria	resident	44	44	1	121	1.86	NA	65	2	100	0.63000	Parafit
Savio_2011	2011	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	horizontal	NA	Erwinia dacicola tephritid symbiont	resident	17	17	10	17	1.00	1.00	17	3	999	0.00700	Parafit
Schardl_et_al_2008	2008	Plant	Plant	Microbe	Fungus	y	Mutualist	Endo	both	NA	Grass endophytic fungi (epichloae)	resident	25	25	16	25	0.96	1.00	26	2	1000	0.00100	Parafit
Sibbald_et_al_2017	2017	Microbe	Amoeba	Microbe	Protist	y	Mutualist	Endo	horizontal	environmental	Paramoeba	resident	7	7	1	7	1.00	1.00	7	1	9999	0.00010	Parafit
Simkova_et_al_2013	2013	Vert	Fish	Invert	Invert	y	Parasite	Ecto	horizontal	autonomous	Monogenea	resident	5	5	1	21	1.00	1.00	21	1	999	0.01300	Parafit
Singh_et_al_2016	2016	Microbe	Fungus	Plant	Plant	y	Mutualist	Endo	horizontal	environmental	Trebouxia algae	resident	23	23	1	28	1.40	1.20	20	1	9999	0.00020	Parafit
Sontowski_et_al_2015	2015	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Wolbachia	resident	14	14	4	14	1.00	1.00	14	1	100000	0.00300	Parafit
Sorenson_et_al_2004	2004	Vert	Bird	Vert	Bird	y	Parasite	Ecto	horizontal	autonomous	Brood parasitic finches	resident	33	33	10	34	1.62	1.43	21	1	1000	0.00100	Parafit
Stireman_et_al_2010	2010	Plant	Plant	Invert	Invert	y	Parasite	Endo	horizontal	autonomous	Gall midge	resident	15	15	10	39	1.39	1.25	28	1	9999	0.00080	Parafit
Sudakaran_et_al_2015	2015	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	NA	NA	Gordonibacter	resident	16	16	7	22	1.00	1.00	22	1	1000	0.97800	Parafit
Sudakaran_et_al_2015	2015	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	NA	NA	Coriobacterium glomerans	resident	14	14	6	26	1.00	1.00	26	1	1000	0.97400	Parafit
Summers_&_Rouse_2014	2014	Invert	Invert	Invert	Invert	y	Parasite	Endo/Ecto	horizontal	autonomous	Myzostomida	resident	53	53	36	78	1.13	1.23	69	10	9999	0.00050	Parafit
Susoy_&_Herrmann_2014	2014	Invert	Invert	Invert	Invert	y	Mutualist	Ecto	vertical	vertical	Nematode	resident	35	35	7	37	1.42	1.31	26	1	9999	0.00010	Parafit
Swafford_&_Bond_2010	2010	Invert	Invert	Invert	Invert	y	Parasite	Ecto	horizontal	contact	Millipede mite	resident	7	7	1	7	1.00	1.00	7	1	9999	0.31900	Parafit
Sweet_et_al_2016	2016	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	Dove wing lice	resident	52	52	25	57	1.33	1.21	43	1	100000	0.00001	Parafit
Sweet_et_al_2016	2016	Vert	Bird	Invert	Invert	y	Parasite	Ecto	both	contact	Dove body lice	resident	52	52	25	58	1.18	1.14	49	4	100000	0.00001	Parafit
Tao_et_al_2013	2013	Vert	Tetrapod	Microbe	Virus	n	Parasite	Endo	horizontal	NA	Polyomavirus	resident	7	7	7	10	1.00	1.00	10	1	99999	0.17592	Parafit
Toju_et_al_2013	2013	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Curculioniphilus	resident	27	27	4	27	1.00	1.00	27	1	99999	0.00001	Parafit
Urban_&_Cryan_2012	2012	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Fulgoroidea associated Betaproteobacteria	resident	30	30	29	30	1.00	1.00	30	1	1000	0.00100	Parafit
Urban_&_Cryan_2012	2012	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Fulgoroidea associated Sulcia	resident	40	40	38	40	1.00	1.00	40	1	1000	0.00300	Parafit
Viale_et_al_2015	2015	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	vertical	vertical	Fruit fly Stammerula	resident	23	23	1	23	1.00	1.00	23	1	9999	0.00100	Parafit
Won_et_al_2008	2008	Invert	Invert	Microbe	Bacterium	n	Mutualist	Endo	horizontal	contact	Mussel thiotroph symbiont	resident	15	15	5	15	1.00	1.00	15	1	100	0.37000	Parafit

A. authors: The authors of the study and the date (citation form).

B. year: Year of publication of the study.

C. host_tax_broad: Separation of the host group according to broader taxonomic units (e.g. vertebrate, invertebrate, microbe, plant).

D. host_tax_fine*: Separation of the host group according to narrower taxonomic units (e.g. fish, tetrapod, bird, invertebrate, protist, bacterium, plant, fungus).

E. symbiont_tax_broad: Separation of the symbiont group according to broader taxonomic units (e.g. vertebrate, invertebrate, microbe, plant).

F. symbiont_tax_fine*: Separation of the symbiont group according to narrower taxonomic units (e.g. invertebrate, protist, virus, bacterium, fungus, plant, bird).

G. symbiont_euk*: Whether the symbiont is eukaryotic (state =‘yes’), or prokaryotic (state=‘no’).

H. symbiosis: The type of symbiont (e.g. parasite or mutualist). For this we followed the definition used by the authors of the study.

I. endo_or_ecto: Whether the symbiont lives outside the host (i.e. is an ectosymbiont), or inside the host (i.e. is an endosymbiont).

J. mode_of_transmission_broad: Whether the symbiont is transmitted vertically, horizontally, or both. For this, we followed the route of transmission specified by the authors of the study

K. mode_of_transmission_fine*: A finer-scale description of the mode of transmission of the symbiont (e.g. contact, vector, bodily fluid, vertical, trophic).

L. symbiont*: Shorthand description of the type of symbiont.

M. Visiting_symbiont?* Whether the symbiont is resident on the host (resident), or makes visits to the host or hosts (visitor).

N. host_tips_linked: The number of individual host taxa included in the cophylogenetic analysis.

O. host_tips_linked_corrected The same measure as for column N, ‘host_tips_linked’, but reduced to only include one member of each host species. This is included because some authors include multiple individuals of the same host species. Without correction, this artificially increases the apparent number of host species included in the study.

P. host_genera: A count of the number of host genera included in the cophylogenetic analysis

Q. total_host_symbioint_links The total number of links between host and symbiont taxa recorded in a study. If all symbionts were strict specialists, this would equal the number of symbionts included in the study. However, because symbionts are often associated with more than one host, this value is often higher than the total number of symbionts included in the study.

R. host_range_link_ratio: An estimation of symbiont host specificity, calculated by dividing the total number of links between hosts and symbionts (i.e. ‘total_host_symbiont_links’, column Q), by the total number of symbionts included in the study (i.e. ‘symbiont_tips_linked’, column T).

S. host_range_taxonomic_breadth: An alternative estimation of symbiont host specificity, calculated by first summing the number of host taxonomic ranks linked to each symbiont (i.e. single host species = 1, multiple host species in the same genera = 2, multiple host genera = 3, multiple host families = 4, multiple host orders = 5), and dividing by the total number of symbionts included in the study (i.e. ‘symbiont_tips_linked’, column T)

T. symbiont_tips_linked The number of individual symbiont taxa included in the cophylogenetic analysis.

U. symbiont_genera: A count of the number of symbiont genera included in the cophylogenetic analysis.

V. no_randomizations: The number of phylogenetic randomizations performed during the cophylogenetic analysis.

W. p_value: The p-value reported for the cophylogenetic analysis, representing the likelihood that host and symbiont phylogenies display cospeciation.

X. method: Whether TreeMap or ParaFit was used to obtain co-divergence or a p value.

Table of sample sizes

We present the number of sample size for two separate methods: TreeMap¹ and ParaFit² (and combined) for effect sizes, papers and different levels of categorical variables (factors).

# selecting out variables, which we used for our analysis
dat <- full_data %>% select(-host_tax_fine, -symbiont_tax_fine, -symbiont_euk, -mode_of_transmission_fine, -symbiont, -`Visiting_symbiont?`)

# making a table of sample sizes for different variables
dat %>% group_by(method) %>% 
  summarise(
    `Effect sizes` = n(),
    Papers = n_distinct(authors),
    `Vertebrate hosts` = sum(host_tax_broad == "Vert", na.rm = T), # na.rm is important when NA exists
    `Invertebrate hosts` = sum(host_tax_broad == "Invert", na.rm = T),
    `Plant hosts`  = sum(host_tax_broad == "Plant", na.rm = T),
    `Microbe hosts` = sum(host_tax_broad == "Microbe", na.rm = T),
    `Vertebrate symbionts` = sum(symbiont_tax_broad  == "Vert", na.rm = T),
    `Invertebrate symbionts` = sum(symbiont_tax_broad  == "Invert", na.rm = T),
    `Plant symbionts`  = sum(symbiont_tax_broad  == "Plant", na.rm = T),
    `Microbe symbionts` = sum(symbiont_tax_broad  == "Microbe", na.rm = T),
    `Parastic relationships` = sum(symbiosis == "Parasite", na.rm = T),
    `Mutualistic relatioships` = sum(symbiosis == "Mutualist", na.rm = T),
    `Ecto-symbionts` = sum(endo_or_ecto  == "Ecto", na.rm = T),
    `Endo-symbionts` = sum(endo_or_ecto == "Endo", na.rm = T),
    `Ecto/endo-symbionts` = sum(endo_or_ecto == "Endo/Ecto", na.rm = T),
    `Horizontal transmission` = sum(mode_of_transmission_broad  == "horizontal", na.rm = T),
    `Vertical transmission` = sum(mode_of_transmission_broad == "vertical", na.rm = T),
    `Horizontal/vertical-transmission` = sum(mode_of_transmission_broad == "both", na.rm = T)
  ) -> n_table1

# transposing the table and creating that table and adding a correct number of the papers for `Combined`
n_authors <- n_distinct(dat$authors) # the total number of papers
n_table2 <-t(n_table1[,-1])
colnames(n_table2) <- n_table1$method
n_table2 %>% as_tibble(rownames = "Number") %>% 
  mutate(Combined = Parafit + TreeMap, Combined = replace(Combined, 2, n_authors)) %>%  
  rename("Number of" = "Number", "ParaFit (n)" = "Parafit", "TreeMap (n)" = "TreeMap",  "Combined (n)" = "Combined") %>% 
  kable() %>% kable_styling("striped", position = "left") %>%
  scroll_box(width = "100%", height = "250px")

Number of	ParaFit (n)	TreeMap (n)	Combined (n)
Effect sizes	140	93	233
Papers	118	78	180
Vertebrate hosts	60	51	111
Invertebrate hosts	39	20	59
Plant hosts	31	18	49
Microbe hosts	10	4	14
Vertebrate symbionts	1	1	2
Invertebrate symbionts	62	49	111
Plant symbionts	3	1	4
Microbe symbionts	74	42	116
Parastic relationships	91	70	161
Mutualistic relatioships	48	23	71
Ecto-symbionts	41	34	75
Endo-symbionts	97	58	155
Ecto/endo-symbionts	2	1	3
Horizontal transmission	84	53	137
Vertical transmission	28	15	43
Horizontal/vertical-transmission	23	25	48

  #pander(split.cell = 40, split.table = Inf) # not as nice as kable

Note that for the numbers of papers do not add up (TreeMap + ParaFit \(\neq\) Combined) because 16 papers used both TreeMap and ParaFit methods (the term “papers” here is our variable auhtors)

Missing data patterns

Below, we have the number of missing data (cells) for all the variables used in the analysis.

# summaring missingness in our dataset
# funs(sum(is.na(.))) needs to be in funs as is.na has "." = each column
dat %>% summarise_all(~sum(is.na(.))) %>% # map(~sum(is.na(.)) # this is an alterantive way 
  t() %>% as_tibble(rownames = "Variable") %>% 
  rename("Number of missing data (n)" = "V1") %>% 
  #pander(split.cell = 40, split.table = Inf)
  kable() %>% kable_styling("striped", position = "left") %>%
  scroll_box(width = "60%", height = "250px")

Variable	Number of missing data (n)
authors	0
year	0
host_tax_broad	0
symbiont_tax_broad	0
symbiosis	1
endo_or_ecto	0
mode_of_transmission_broad	5
host_tips_linked	0
host_tips_linked_corrected	0
host_genera	6
total_host_symbioint_links	3
host_range_link_ratio	3
host_range_taxonomic_breadth	7
symbiont_tips_linked	0
symbiont_genera	9
no_randomizations	0
p_value	0
method	0

# an alternative method using the mi package
#missing_data_tbl <- missing_data.frame(as.data.frame(data))
#show(missing_data_tbl) 

Meta-analysis

Calculating effect sizes

Here, we created our effect size (correlation coefficient r and its Fisher’s z transformation Zr) from p values and associated sample sizes³. We used the sum of host_tips_linked_corrected and symbiont_tips_linked as our sample size (i.e., the number of both host and symbiont species) for each effect size (an indicator of congruence). Also, we created a column with unique ID for each observation (i.e. an observation level random effect), termed observation, which is required for the rma.mv function in metafor⁴.

dat %<>% # getting sample size & observation level random effect
mutate(., sample_size = host_tips_linked_corrected + symbiont_tips_linked, observation = factor(1:nrow(.)))
# calcuating effect size somehow it did not run by piping with above (not
# recognising the first aurgument)
dat %<>% p_to_Zr(p_value, sample_size)

Meta-analytic model: testing Fahrenholz’s rule

First, we checked what random effects should be put into the main model to do this we fitted two random effects, authors (i.e. study IDs) and observation; the former term was added to account for non-independence of effect sizes originating from the same papers (i.e., authors).

# 2 random effects & model AIC note that probably only base stuff works
# outside of main chunck so need to create AIC here
ma_test1 <- rma.mv(yi = Zr, V = VZr, random = list(~1 | authors, ~1 | observation), 
    data = dat)
aic1 <- AIC(ma_test1)

# 1 random effect & model AIC
ma_test2 <- rma.mv(yi = Zr, V = VZr, random = ~1 | authors, data = dat)
aic2 <- AIC(ma_test2)

The model (ma_test1), which included both random factors, had a larger AIC value (39.91) than the model with only one random effect (37.91) . This is because observation hardly accounted any variance (< 0.0001) compared to authors (0.0327). Therefore, we only had authors as our random factor in subsequent analyses.

We ran intercept models (meta-analyses) with 3 different datasets (ParaFit, TreeMap and both combined; see the explanation of method above). Also, we note that we used adjustments for test statistics and confidence intervals (test = "t"), which is similar to (but not the same as) those proposed by Kanpp and Hartung⁵; probably this approach is a more conservative approach.

Running Multilevel Meta-analytic models with 3 datasets

# think about making this into a tibble meta-analysis with Parafit
ma_parafit <- rma.mv(yi = Zr, V = VZr, random = ~1 | authors, test = "t", subset = which(method == 
    "TreeMap"), data = dat)

# meta-analysis with TreeMap
ma_treemap <- rma.mv(yi = Zr, V = VZr, random = ~1 | authors, test = "t", subset = which(method == 
    "Parafit"), data = dat)

# meta-analysis with all the data combined
ma_all <- rma.mv(yi = Zr, V = VZr, test = "t", random = ~1 | authors, data = dat)

Supplementary Table 1: Overall effects (meta-analytic means), 95% confidence intervals (CIs), variance components (V) and heterogeneity, I² (I2)⁶ from metafor model using the 3 datasets (ParaFit, TreeMap and both combined, or All). Note that in these models, I²_[total] = I²_[authors] (see^7,8), as we only have one random factor.

# getting I2 for the models could use map()
i2_treemap <- I2(ma_treemap)
i2_parafit <- I2(ma_parafit)
i2_all <- I2(ma_all)
# creating a table
tibble(Dataset = c("Parafit", "TreeMap", "All"), `Overall mean (Zr)` = c(ma_parafit$b, 
    ma_treemap$b, ma_all$b), `Lower CI [0.025]` = c(ma_parafit$ci.lb, ma_treemap$ci.lb, 
    ma_all$ci.lb), `Upper C  [0.975]` = c(ma_parafit$ci.ub, ma_treemap$ci.ub, 
    ma_all$ci.ub), `V[authors]` = c(ma_parafit$sigma2, ma_treemap$sigma2, ma_all$sigma2), 
    `I2[total]` = c(i2_parafit[1], i2_treemap[1], i2_all[1])) %>% kable("html", 
    digits = 3) %>% kable_styling("striped", position = "left")

Dataset	Overall mean (Zr)	Lower CI [0.025]	Upper C [0.975]	V[authors]	I2[total]
Parafit	0.333	0.273	0.394	0.036	0.523
TreeMap	0.350	0.308	0.392	0.028	0.600
All	0.345	0.309	0.381	0.033	0.590

These models gave all consistent results including heterogeneity. Given these results, we proceeded with only analyzing the whole dataset (All) from this on.

# https://stackoverflow.com/questions/41919023/ggplot-adding-image-on-top-right-in-two-plots-with-different-scales
# how to add png files to the figure (above) reading image
image_mutualism <- readPNG("../images/mutualism_transparentbg.png")
image_parasitism <- readPNG("../images/parasitism_transparentbg.png")

# creating a table of results
pred_ma <- get_pred(ma_all)
effect_ma <- get_est(ma_all) %>% left_join(pred_ma)

# creating a forest plot
fig_ma <- ggplot(data = effect_ma, aes(x = tanh(estimate), y = "Overall mean")) + 
    scale_x_continuous(limits = c(-1, 1), breaks = seq(-1, 1, by = 0.2)) + geom_quasirandom(data = dat, 
    aes(x = tanh(Zr), y = "Overall mean", size = (1/VZr) + 3), groupOnX = FALSE, 
    alpha = 0.2) + # precition interval (PI)
geom_errorbarh(aes(xmin = tanh(lowerPR), xmax = tanh(upperPR)), height = 0, 
    show.legend = F, size = 0.5, alpha = 0.6) + # CI
geom_errorbarh(aes(xmin = tanh(lowerCL), xmax = tanh(upperCL)), height = 0, 
    show.legend = F, size = 1.2) + 
geom_vline(xintercept = 0, linetype = 2, colour = "black", alpha = 0.3) + # creating dots and different size (bee-swarm and bubbles)
geom_point(size = 3, shape = 21, fill = "black") + annotate("text", x = 0.93, 
    y = 1.15, label = paste("italic(k)==", length(dat$Zr)), parse = TRUE, hjust = "left", 
    size = 3.5) + labs(x = expression(paste(italic(r), " (correlation)")), y = "", 
    size = expression(paste(italic(n), " (# of species)"))) + theme_bw() + theme(legend.position = c(0, 
    1), legend.justification = c(0, 1)) + theme(legend.direction = "horizontal") + 
    # theme(legend.background = element_rect(fill = 'white', colour = 'black'))
# +
theme(legend.background = element_blank()) + theme(axis.text.y = element_text(size = 10, 
    colour = "black", hjust = 0.5, angle = 90)) + annotation_custom(rasterGrob(image_mutualism), 
    xmin = -1.1, xmax = -0.9, ymin = 0.6, ymax = 1.2) + annotation_custom(rasterGrob(image_parasitism), 
    xmin = -0.9, xmax = -0.7, ymin = 0.6, ymax = 1.2)
# ggsave(plot = fig_ma, filename = 'fig_2a.pdf', height = 2, width = 8)
# ggploty 0 does not work (Error in unique.default(x) : unimplemented type
# 'expression' in 'HashTableSetup')
fig_ma

# for Fig 3

a <- ggplot(data = effect_ma, aes(x = tanh(estimate), y = "Overall mean")) + 
    scale_x_continuous(limits = c(-1, 1), breaks = seq(-1, 1, by = 0.2)) + geom_quasirandom(data = dat, 
    aes(x = tanh(Zr), y = "Overall mean", size = (1/VZr) + 3), groupOnX = FALSE, 
    alpha = 0.2) + # precition interval (PI)
geom_errorbarh(aes(xmin = tanh(lowerPR), xmax = tanh(upperPR)), height = 0, 
    show.legend = F, size = 0.5, alpha = 0.6) + # CI
geom_errorbarh(aes(xmin = tanh(lowerCL), xmax = tanh(upperCL)), height = 0, 
    show.legend = F, size = 1.2) + 
geom_vline(xintercept = 0, linetype = 2, colour = "black", alpha = 0.3) + # creating dots and different size (bee-swarm and bubbles)
geom_point(size = 3, shape = 21, fill = "black") + annotate("text", x = 0.93, 
    y = 1.15, label = paste("italic(k)==", length(dat$Zr)), parse = TRUE, hjust = "left", 
    size = 3.5) + labs(x = "", y = "", size = expression(paste(italic(n), " (# of species)")), 
    tag = "a") + theme_bw() + theme(legend.position = c(0, 1), legend.justification = c(0, 
    1)) + theme(legend.direction = "horizontal") + # theme(legend.background = element_rect(fill = 'white', colour = 'black'))
# +
theme(legend.background = element_blank()) + theme(axis.text.y = element_text(size = 10, 
    colour = "black", hjust = 0.5, angle = 90)) + annotation_custom(rasterGrob(image_mutualism), 
    xmin = -1.1, xmax = -0.9, ymin = 0.6, ymax = 1.2) + annotation_custom(rasterGrob(image_parasitism), 
    xmin = -0.9, xmax = -0.7, ymin = 0.6, ymax = 1.2)

Figure 3a: A forest plot showing the meta-analytic mean (mean effect size) with its 95% confidence interval (thick line) and 95% prediction interval (thin line) with observed effect sizes based on various sample sizes.

Meta-regression

We ran a univariate meta-regression model for each of following moderators: 1) symbiosis, 2) host_tax_broad, 3) symbiont_tax_broad, 4) host_range_link_ratio, 5) host_range_taxonomic_breadth, 6) mode_of_transmission_broad, and 7) endo_or_ecto. The results from these models are presented in the main text.

In addition to these, we ran three more univariate models: 1) host_tax_symbiosis (equivalent to the interaction term between symbiosis and host_tax_symbiosis; symbiosis*host_tax_symbiosis), 2) symbiont_tax_symbiosis (symbiosis*symbiont_tax_broad), 3) host_symbiont_tax (host_tax_symbiosis*symbiont_tax_broad) and 4) symbiosis_transmission (symbiosis*mode_of_transmission_broad). There moderators are created below:

dat %<>% 
    # host_tax_broad*symbiosis (host_tax_symbiosis) 
  mutate(host_tax_symbiosis = str_c(host_tax_broad, symbiosis), 
         host_tax_symbiosis = ifelse(host_tax_symbiosis == "InvertNA", NA, host_tax_symbiosis),
         host_tax_symbiosis = factor(host_tax_symbiosis),
         # symbiont_tax_broad*symbiosis (symbiont_tax_symbiosis)     
         symbiont_tax_symbiosis = factor(str_c(symbiont_tax_broad, symbiosis)),
         # host_tax_broad*symbiont_tax_broad (host_symbiont_tax)     
         host_symbiont_tax  = factor(str_c(host_tax_broad, symbiont_tax_broad)),
         # symbiosis*mode_of_transmission_broad (symbiosis_transmission)
         symbiosis_transmission  = factor(str_c(symbiosis, mode_of_transmission_broad)),
         # whether p values were the smallest value given the number of randamization - limit_researched (Yes = 1, No = 0)
         limit_rearched = if_else(abs((1/p_value) - no_randomizations) <= 1, 1, 0))

Univariate (uni-predictor) analyses

We first conducted a series of meta-regression model with one predictor.

The type of symbiosis: parasitism vs. mutualism

# meta-regression: mutiple intercepts
mr_symbiosis1 <- rma.mv(yi = Zr, V = VZr, mods = ~symbiosis - 1, test = "t", 
    random = ~1 | authors, data = dat)
# meta-regression: contrast
mr_symbiosis2 <- rma.mv(yi = Zr, V = VZr, mods = ~symbiosis, test = "t", random = ~1 | 
    authors, data = dat)

Supplementary Table 2: Regression coefficients (Estimate), 95% confidence intervals (CIs), variance components (V) and variance explained, R²_[marginal]⁹ (R2) from the meta-regression with symbiosis.

# getting marginal R2
r2_symbiosis1 <- R2(mr_symbiosis1)

# getting estimates
res_symbiosis1 <- get_est(mr_symbiosis1, mod = "symbiosis")
res_symbiosis2 <- get_est(mr_symbiosis2, mod = "symbiosis")

# creating a table
tibble(`Fixed effect` = c(as.character(res_symbiosis1$name), cont_gen(res_symbiosis1$name)), 
    Estimate = c(res_symbiosis1$estimate, res_symbiosis2$estimate[2]), `Lower CI [0.025]` = c(res_symbiosis1$lowerCL, 
        res_symbiosis2$lowerCL[2]), `Upper CI  [0.975]` = c(res_symbiosis1$upperCL, 
        res_symbiosis2$upperCL[2]), `V[authors]` = c(mr_symbiosis1$sigma2, rep(NA, 
        2)), R2 = c(r2_symbiosis1[1], rep(NA, 2))) %>% kable("html", digits = 3) %>% 
    kable_styling("striped", position = "left")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	V[authors]	R2
Mutualist	0.411	0.349	0.473	0.031	0.064
Parasite	0.311	0.268	0.354	NA	NA
Mutualist-Parasite	-0.100	-0.174	-0.026	NA	NA

# adding sample size (k) for each category
k_symbiosis <- dat %>% group_by(symbiosis) %>% count()
# getting estimates and predicitons
pred_symbiosis <- get_pred(mr_symbiosis1, mod = "symbiosis") 
res_symbiosis1 <- left_join(res_symbiosis1, k_symbiosis, by =  c("name" = "symbiosis"))  %>% left_join(pred_symbiosis)
#res_symbiosis1 
# drawing a funnel plot - fig 2b
fig_symbiosis <- ggplot(data = res_symbiosis1, aes(x = tanh(estimate), y = name)) +
  scale_x_continuous(limits=c(-1, 1), breaks = seq(-1, 1, by = 0.2) ) +
  geom_quasirandom(data = dat %>% filter(!is.na(symbiosis)), 
                   aes(x= tanh(Zr), y = symbiosis, size = ((1/VZr) + 3), colour = symbiosis), groupOnX = FALSE, alpha=0.2) + 
  # 95 %precition interval (PI)
  geom_errorbarh(aes(xmin = tanh(lowerPR), xmax = tanh(upperPR)),  height = 0, show.legend = F, size = 0.5, alpha = 0.6) +
  # 95 %CI
  geom_errorbarh(aes(xmin = tanh(lowerCL), xmax = tanh(upperCL)),  height = 0, show.legend = F, size = 1.2) +
  geom_vline(xintercept = 0, linetype = 2, colour = "black", alpha = 0.3) +
  # creating dots and different size (bee-swarm and bubbles)
  geom_point(aes(fill = name), size = 3, shape = 21) + #
  # setting colours
  scale_color_manual(values = c("Mutualist" = "#E69F00", "Parasite" = "#56B4E9")) +
  scale_fill_manual(values = c("Mutualist" = "#E69F00", "Parasite" = "#56B4E9")) +
  annotate('text', x = 0.93, y = c(1.15, 2.15), label= paste("italic(k)==", res_symbiosis1$n), parse = TRUE, hjust = "left", size = 3.5) +
  labs(x = expression(paste(italic(r), " (correlation)")), y = "", size = expression(paste(italic(n), " (# of species)")) ) +
  guides(fill = "none", colour = "none") +
  theme_bw() +
  theme(legend.position= c(0, 1), legend.justification = c(0,1)) +
  theme(legend.direction="horizontal") +
  #theme(legend.background = element_rect(fill = "white", colour = "black")) +
  theme(legend.background = element_blank()) +
  theme(axis.text.y = element_text(size = 10, colour ="black", hjust = 0.5, angle = 90)) +
  # putting pictures in
  annotation_custom(rasterGrob(image_mutualism), xmin = -1, xmax = -0.8, ymin = 0.6, ymax = 1.2) + 
  annotation_custom(rasterGrob(image_parasitism), xmin = -1, xmax = -0.8, ymin = 1.6, ymax = 2.2)

fig_symbiosis

# fig 3

b <- ggplot(data = res_symbiosis1, aes(x = tanh(estimate), y = name)) +
  scale_x_continuous(limits=c(-1, 1), breaks = seq(-1, 1, by = 0.2) ) +
  geom_quasirandom(data = dat %>% filter(!is.na(symbiosis)), 
                   aes(x= tanh(Zr), y = symbiosis, size = ((1/VZr) + 3), colour = symbiosis), groupOnX = FALSE, alpha=0.2) + 
  # 95 %precition interval (PI)
  geom_errorbarh(aes(xmin = tanh(lowerPR), xmax = tanh(upperPR)),  height = 0, show.legend = F, size = 0.5, alpha = 0.6) +
  # 95 %CI
  geom_errorbarh(aes(xmin = tanh(lowerCL), xmax = tanh(upperCL)),  height = 0, show.legend = F, size = 1.2) +
  geom_vline(xintercept = 0, linetype = 2, colour = "black", alpha = 0.3) +
  # creating dots and different size (bee-swarm and bubbles)
  geom_point(aes(fill = name), size = 3, shape = 21) + #
  # setting colours
  scale_color_manual(values = c("Mutualist" = "#E69F00", "Parasite" = "#56B4E9")) +
  scale_fill_manual(values = c("Mutualist" = "#E69F00", "Parasite" = "#56B4E9")) +
  annotate('text', x = 0.93, y = c(1.15, 2.15), label= paste("italic(k)==", res_symbiosis1$n), parse = TRUE, hjust = "left", size = 3.5) +
  labs(x = "", y = "", tag = "b") +
  guides(fill = "none", colour = "none") +
  theme_bw() +
  theme(legend.position="none") +
  theme(axis.text.y = element_text(size = 10, colour ="black",hjust = 0.5, angle = 90)) +
  # putting pictures in
  annotation_custom(rasterGrob(image_mutualism), xmin = -1, xmax = -0.8, ymin = 0.6, ymax = 1.2) + 
  annotation_custom(rasterGrob(image_parasitism), xmin = -1, xmax = -0.8, ymin = 1.6, ymax = 2.2)

Figure 3b: A forest plot showing the group-wise means (the categorical variable symbiosis) with their 95% confidences interval (thick lines) and 95% prediction intervals (thin lines) with observed effect sizes based on various sample sizes.

The effect of host taxa

# reordering
dat$host_tax_broad <- factor(dat$host_tax_broad, levels = c("Microbe", "Plant", 
    "Invert", "Vert"))

# meta-regression: mutiple intercepts
mr_host_tax_broad1 <- rma.mv(yi = Zr, V = VZr, mods = ~host_tax_broad - 1, test = "t", 
    random = ~1 | authors, data = dat)

# meta-regression: contrast 1
mr_host_tax_broad2 <- rma.mv(yi = Zr, V = VZr, mods = ~host_tax_broad, test = "t", 
    random = ~1 | authors, data = dat)

# meta-regression: contrast 2
mr_host_tax_broad3 <- rma.mv(yi = Zr, V = VZr, mods = ~relevel(host_tax_broad, 
    ref = "Plant"), test = "t", random = ~1 | authors, data = dat)

# meta-regression: contrast 3
mr_host_tax_broad4 <- rma.mv(yi = Zr, V = VZr, mods = ~relevel(host_tax_broad, 
    ref = "Invert"), test = "t", random = ~1 | authors, data = dat)

Supplementary Table 3: Regression coefficients (Estimate), 95% confidence intervals (CIs), variance components (V) and variance explained, R²_[marginal] (R2) from the meta-regression with host_tax_broad.

# getting marginal R2
r2_host_tax_broad1 <- R2(mr_host_tax_broad1)

# getting estimates
res_host_tax_broad1 <- get_est(mr_host_tax_broad1, mod = "host_tax_broad")
res_host_tax_broad2 <- get_est(mr_host_tax_broad2, mod = "host_tax_broad")
# the name bit does not work if relevel....
res_host_tax_broad3 <- get_est(mr_host_tax_broad3, mod = "host_tax_broad")
res_host_tax_broad4 <- get_est(mr_host_tax_broad4, mod = "host_tax_broad")

# creating a table
tibble(`Fixed effect` = c(as.character(res_host_tax_broad1$name), cont_gen(res_host_tax_broad1$name)), 
    Estimate = c(res_host_tax_broad1$estimate, res_host_tax_broad2$estimate[-1], 
        res_host_tax_broad3$estimate[-(1:2)], res_host_tax_broad4$estimate[-(1:3)]), 
    `Lower CI [0.025]` = c(res_host_tax_broad1$lowerCL, res_host_tax_broad2$lowerCL[-1], 
        res_host_tax_broad3$lowerCL[-(1:2)], res_host_tax_broad4$lowerCL[-(1:3)]), 
    `Upper CI  [0.975]` = c(res_host_tax_broad1$upperCL, res_host_tax_broad2$upperCL[-1], 
        res_host_tax_broad3$upperCL[-(1:2)], res_host_tax_broad4$upperCL[-(1:3)]), 
    `V[authors]` = c(mr_host_tax_broad1$sigma2, rep(NA, 9)), R2 = c(r2_host_tax_broad1[1], 
        rep(NA, 9))) %>% kable("html", digits = 3) %>% kable_styling("striped", 
    position = "left") %>% scroll_box(width = "100%", height = "300px")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	V[authors]	R2
Microbe	0.557	0.404	0.710	0.03	0.168
Plant	0.244	0.168	0.321	NA	NA
Invert	0.412	0.344	0.480	NA	NA
Vert	0.327	0.277	0.378	NA	NA
Microbe-Plant	-0.312	-0.483	-0.141	NA	NA
Microbe-Invert	-0.144	-0.312	0.023	NA	NA
Microbe-Vert	-0.229	-0.390	-0.068	NA	NA
Plant-Invert	0.168	0.065	0.270	NA	NA
Plant-Vert	0.083	-0.009	0.175	NA	NA
Invert-Vert	-0.085	-0.168	-0.002	NA	NA

# getting images
image_invertebrate_host <- readPNG("../images/invertebrate_host_transparentbg.png")
image_microbe_host <- readPNG("../images/microbe_host_transparentbg.png")
image_vertebrate_host <- readPNG("../images/vertebrate_host_transparentbg.png")
image_plant_host <- readPNG("../images/plant_host_transparentbg.png")

# adding sample size (k) for each category
k_host_tax_broad <- dat %>% group_by(host_tax_broad) %>% count()
# getting estimates and predicitons
pred_host_tax_broad <- get_pred(mr_host_tax_broad1, mod = "host_tax_broad") 
res_host_tax_broad1 <- left_join(res_host_tax_broad1, k_host_tax_broad, by =  c("name" = "host_tax_broad"))  %>% left_join(pred_host_tax_broad)
#res_symbiosis1 
# drawing a funnel plot - fig 2b
fig_host_tax_broad <- ggplot(data = res_host_tax_broad1, aes(x = tanh(estimate), y = name)) +
  scale_x_continuous(limits=c(-1, 1), breaks = seq(-1, 1, by = 0.2) ) +
  geom_quasirandom(data = dat %>% filter(!is.na(host_tax_broad)), 
                   aes(x= tanh(Zr), y = host_tax_broad, size = ((1/VZr) + 3), colour = host_tax_broad), groupOnX = FALSE, alpha=0.4) + 
  # 95 %precition interval (PI)
  geom_errorbarh(aes(xmin = tanh(lowerPR), xmax = tanh(upperPR)),  height = 0, show.legend = F, size = 0.5, alpha = 0.6) +
  # 95 %CI
  geom_errorbarh(aes(xmin = tanh(lowerCL), xmax = tanh(upperCL)),  height = 0, show.legend = F, size = 1.2) +
  geom_vline(xintercept = 0, linetype = 2, colour = "black", alpha = 0.3) +
  # creating dots and different size (bee-swarm and bubbles)
  geom_point(aes(fill = name), size = 3, shape = 21) + #
  # setting colours
  scale_color_manual(values = c("Microbe" = "#009E73",  "Plant" = "#F0E422",  "Invert"= "#0072B2",  "Vert" = "#D55E00")) +
  scale_fill_manual(values = c("Microbe" = "#009E73",  "Plant" = "#F0E422",  "Invert"= "#0072B2",  "Vert" = "#D55E00")) +
  scale_y_discrete(labels = c("Microbe" = "Microbe",  "Plant" = "Plant",  "Invert"= "Invertebrate",  "Vert" = "Vertebrate")) +
  annotate('text', x = 0.93, y = 1:4 + 0.15, label= paste("italic(k)==", res_host_tax_broad1$n), parse=TRUE, hjust = "left", size=3.5) +
  labs(x = expression(paste(italic(r), " (correlation)")), y = "", size = expression(paste(italic(n), " (# of species)")) ) +
  guides(fill = "none", colour = "none") +
  theme_bw() +
  theme(legend.position= c(0, 1), legend.justification = c(0,1)) +
  theme(legend.direction="horizontal") +
  #theme(legend.background = element_rect(fill = "white", colour = "black")) +
  theme(legend.background = element_blank()) +
  theme(axis.text.y = element_text(size = 10, colour ="black", hjust = 0.5, angle = 90)) +
  # putting pictures in
  annotation_custom(rasterGrob(image_microbe_host), xmin = -1, xmax = -0.8, ymin = 0.6, ymax = 1.2) + 
  annotation_custom(rasterGrob(image_plant_host), xmin = -1, xmax = -0.8, ymin = 1.6, ymax = 2.2) +
  annotation_custom(rasterGrob(image_invertebrate_host), xmin = -1, xmax = -0.8, ymin = 2.6, ymax = 3.2) + 
  annotation_custom(rasterGrob(image_vertebrate_host), xmin = -1, xmax = -0.8, ymin = 3.6, ymax = 4.2)

fig_host_tax_broad

# fig 3c
c <- ggplot(data = res_host_tax_broad1, aes(x = tanh(estimate), y = name)) +
  scale_x_continuous(limits=c(-1, 1), breaks = seq(-1, 1, by = 0.2) ) +
  geom_quasirandom(data = dat %>% filter(!is.na(host_tax_broad)), 
                   aes(x= tanh(Zr), y = host_tax_broad, size = ((1/VZr) + 3), colour = host_tax_broad), groupOnX = FALSE, alpha=0.4) + 
  # 95 %precition interval (PI)
  geom_errorbarh(aes(xmin = tanh(lowerPR), xmax = tanh(upperPR)),  height = 0, show.legend = F, size = 0.5, alpha = 0.6) +
  # 95 %CI
  geom_errorbarh(aes(xmin = tanh(lowerCL), xmax = tanh(upperCL)),  height = 0, show.legend = F, size = 1.2) +
  geom_vline(xintercept = 0, linetype = 2, colour = "black", alpha = 0.3) +
  # creating dots and different size (bee-swarm and bubbles)
  geom_point(aes(fill = name), size = 3, shape = 21) + #
  # setting colours
  scale_color_manual(values = c("Microbe" = "#009E73",  "Plant" = "#F0E422",  "Invert"= "#0072B2",  "Vert" = "#D55E00")) +
  scale_fill_manual(values = c("Microbe" = "#009E73",  "Plant" = "#F0E422",  "Invert"= "#0072B2",  "Vert" = "#D55E00")) +
  scale_y_discrete(labels = c("Microbe" = "Microbe",  "Plant" = "Plant",  "Invert"= "Invertebrate",  "Vert" = "Vertebrate")) +
  annotate('text', x = 0.93, y = 1:4 + 0.15, label= paste("italic(k)==", res_host_tax_broad1$n), parse=TRUE, hjust = "left", size=3.5) +
  labs(x = "", y = "", size = expression(paste(italic(n), " (# of species)")) , tag = "c") +
  guides(fill = "none", colour = "none") +
  theme_bw() +
  theme(legend.position="none") +
  theme(axis.text.y = element_text(size = 10, colour ="black",hjust = 0.5, angle = 90)) +
  # putting pictures in
  annotation_custom(rasterGrob(image_microbe_host), xmin = -1, xmax = -0.8, ymin = 0.6, ymax = 1.2) + 
  annotation_custom(rasterGrob(image_plant_host), xmin = -1, xmax = -0.8, ymin = 1.6, ymax = 2.2) +
  annotation_custom(rasterGrob(image_invertebrate_host), xmin = -1, xmax = -0.8, ymin = 2.6, ymax = 3.2) + 
  annotation_custom(rasterGrob(image_vertebrate_host), xmin = -1, xmax = -0.8, ymin = 3.6, ymax = 4.2)

Figure 3c: A forest plot showing the group-wise means (the categorical variable host_tax_broad) with their 95% confidences interval (thick lines) and 95% prediction intervals (thin lines) with observed effect sizes based on various sample sizes.

The effect of symbiont taxa

# reordering
dat$symbiont_tax_broad <- factor(dat$symbiont_tax_broad, levels = c("Microbe", 
    "Plant", "Invert", "Vert"))

# sizes <- factor(sizes, levels = c('small', 'medium', 'large')) sizes > [1]
# small large large small medium > Levels: small medium large
# meta-regression: mutiple intercepts
mr_symbiont_tax_broad1 <- rma.mv(yi = Zr, V = VZr, mods = ~symbiont_tax_broad - 
    1, test = "t", random = ~1 | authors, data = dat)

# meta-regression: contrast 1
mr_symbiont_tax_broad2 <- rma.mv(yi = Zr, V = VZr, mods = ~symbiont_tax_broad, 
    test = "t", random = ~1 | authors, data = dat)

# meta-regression: contrast 2
mr_symbiont_tax_broad3 <- rma.mv(yi = Zr, V = VZr, mods = ~relevel(symbiont_tax_broad, 
    ref = "Plant"), test = "t", random = ~1 | authors, data = dat)

# meta-regression: contrast 3
mr_symbiont_tax_broad4 <- rma.mv(yi = Zr, V = VZr, mods = ~relevel(symbiont_tax_broad, 
    ref = "Invert"), test = "t", random = ~1 | authors, data = dat)

Supplementary Table 4: Regression coefficients (Estimate), 95% confidence intervals (CIs), variance components (V) and variance explained, R²_[marginal] (R2) from the meta-regression with symbiont_tax_broad.

# getting marginal R2
r2_symbiont_tax_broad1 <- R2(mr_symbiont_tax_broad1)

# getting estimates
res_symbiont_tax_broad1 <- get_est(mr_symbiont_tax_broad1, mod = "symbiont_tax_broad")
res_symbiont_tax_broad2 <- get_est(mr_symbiont_tax_broad2, mod = "symbiont_tax_broad")
res_symbiont_tax_broad3 <- get_est(mr_symbiont_tax_broad3, mod = "symbiont_tax_broad")
res_symbiont_tax_broad4 <- get_est(mr_symbiont_tax_broad4, mod = "symbiont_tax_broad")

# creating a table
tibble(`Fixed effect` = c(as.character(res_symbiont_tax_broad1$name), cont_gen(res_symbiont_tax_broad1$name)), 
    Estimate = c(res_symbiont_tax_broad1$estimate, res_symbiont_tax_broad2$estimate[-1], 
        res_symbiont_tax_broad3$estimate[-(1:2)], res_symbiont_tax_broad4$estimate[-(1:3)]), 
    `Lower CI [0.025]` = c(res_symbiont_tax_broad1$lowerCL, res_symbiont_tax_broad2$lowerCL[-1], 
        res_symbiont_tax_broad3$lowerCL[-(1:2)], res_symbiont_tax_broad4$lowerCL[-(1:3)]), 
    `Upper CI  [0.975]` = c(res_symbiont_tax_broad1$upperCL, res_symbiont_tax_broad2$upperCL[-1], 
        res_symbiont_tax_broad3$upperCL[-(1:2)], res_symbiont_tax_broad4$upperCL[-(1:3)]), 
    `V[authors]` = c(mr_symbiont_tax_broad1$sigma2, rep(NA, 9)), R2 = c(r2_symbiont_tax_broad1[1], 
        rep(NA, 9))) %>% kable("html", digits = 3) %>% kable_styling("striped", 
    position = "left") %>% scroll_box(width = "100%", height = "300px")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	V[authors]	R2
Microbe	0.353	0.303	0.403	0.032	0.082
Plant	0.733	0.440	1.027	NA	NA
Invert	0.324	0.270	0.377	NA	NA
Vert	0.291	-0.112	0.694	NA	NA
Microbe-Plant	0.380	0.083	0.678	NA	NA
Microbe-Invert	-0.029	-0.102	0.044	NA	NA
Microbe-Vert	-0.062	-0.468	0.345	NA	NA
Plant-Invert	-0.410	-0.708	-0.111	NA	NA
Plant-Vert	-0.442	-0.941	0.057	NA	NA
Invert-Vert	-0.033	-0.440	0.374	NA	NA

# getting images
image_invertebrate_parasite <- readPNG("../images/invertebrate_parasite_transparentbg.png")
image_microbe_parasite <- readPNG("../images/microbe_parasite_transparentbg.png")
image_vertebrate_parasite <- readPNG("../images/vertebrate_parasite_transparentbg.png")
image_plant_parasite <- readPNG("../images/plant_parasite_transparentbg.png")

# adding sample size (k) for each category
k_symbiont_tax_broad <- dat %>% group_by(symbiont_tax_broad) %>% count()
# getting estimates and predicitons
pred_symbiont_tax_broad <- get_pred(mr_symbiont_tax_broad1, mod = "symbiont_tax_broad") 
res_symbiont_tax_broad1 <- left_join(res_symbiont_tax_broad1, k_symbiont_tax_broad, by =  c("name" = "symbiont_tax_broad"))  %>% left_join(pred_symbiont_tax_broad)
#res_symbiosis1 
# drawing a funnel plot - fig 2b
fig_symbiont_tax_broad <- ggplot(data = res_symbiont_tax_broad1, aes(x = tanh(estimate), y = name)) +
  scale_x_continuous(limits=c(-1, 1), breaks = seq(-1, 1, by = 0.2) ) +
  geom_quasirandom(data = dat %>% filter(!is.na(symbiont_tax_broad)), 
                   aes(x= tanh(Zr), y = symbiont_tax_broad, size = ((1/VZr) + 3), colour = symbiont_tax_broad), groupOnX = FALSE, alpha=0.4) + 
  # 95 %precition interval (PI)
  geom_errorbarh(aes(xmin = tanh(lowerPR), xmax = tanh(upperPR)),  height = 0, show.legend = F, size = 0.5, alpha = 0.6) +
  # 95 %CI
  geom_errorbarh(aes(xmin = tanh(lowerCL), xmax = tanh(upperCL)),  height = 0, show.legend = F, size = 1.2) +
  geom_vline(xintercept = 0, linetype = 2, colour = "black", alpha = 0.3) +
  # creating dots and different size (bee-swarm and bubbles)
  geom_point(aes(fill = name), size = 3, shape = 21) + #
  # setting colours
  scale_color_manual(values = c("Microbe" = "#009E73",  "Plant" = "#F0E422",  "Invert"= "#0072B2",  "Vert" = "#D55E00" )) +
  scale_fill_manual(values = c("Microbe" = "#009E73",  "Plant" = "#F0E422",  "Invert"= "#0072B2",  "Vert" = "#D55E00")) +
  scale_y_discrete(labels = c("Microbe" = "Microbe",  "Plant" = "Plant",  "Invert"= "Invertebrate",  "Vert" = "Vertebrate")) +
  annotate('text', x = 0.93, y = 1:4 + 0.15, label= paste("italic(k)==", res_symbiont_tax_broad1$n), parse=TRUE, hjust = "left", size=3.5) +
  labs(x = expression(paste(italic(r), " (correlation)")), y = "", size = expression(paste(italic(n), " (# of species)")) ) +
  guides(fill = "none", colour = "none") +
  theme_bw() +
  theme(legend.position= c(0, 1), legend.justification = c(0,1)) +
  theme(legend.direction="horizontal") +
  #theme(legend.background = element_rect(fill = "white", colour = "black")) +
  theme(legend.background = element_blank()) +
  theme(axis.text.y = element_text(size = 10, colour ="black", hjust = 0.5, angle = 90)) +
  # putting pictures in
  annotation_custom(rasterGrob(image_microbe_parasite), xmin = -1, xmax = -0.8, ymin = 0.6, ymax = 1.2) + 
  annotation_custom(rasterGrob(image_plant_parasite), xmin = -1, xmax = -0.8, ymin = 1.6, ymax = 2.2) +
  annotation_custom(rasterGrob(image_invertebrate_parasite), xmin = -1, xmax = -0.8, ymin = 2.6, ymax = 3.2) + 
  annotation_custom(rasterGrob(image_vertebrate_parasite), xmin = -1, xmax = -0.8, ymin = 3.6, ymax = 4.2)

fig_symbiont_tax_broad

# fig 3d
d <- ggplot(data = res_symbiont_tax_broad1, aes(x = tanh(estimate), y = name)) +
  scale_x_continuous(limits=c(-1, 1), breaks = seq(-1, 1, by = 0.2) ) +
  geom_quasirandom(data = dat %>% filter(!is.na(symbiont_tax_broad)), 
                   aes(x= tanh(Zr), y = symbiont_tax_broad, size = ((1/VZr) + 3), colour = symbiont_tax_broad), groupOnX = FALSE, alpha=0.4) + 
  # 95 %precition interval (PI)
  geom_errorbarh(aes(xmin = tanh(lowerPR), xmax = tanh(upperPR)),  height = 0, show.legend = F, size = 0.5, alpha = 0.6) +
  # 95 %CI
  geom_errorbarh(aes(xmin = tanh(lowerCL), xmax = tanh(upperCL)),  height = 0, show.legend = F, size = 1.2) +
  geom_vline(xintercept = 0, linetype = 2, colour = "black", alpha = 0.3) +
  # creating dots and different size (bee-swarm and bubbles)
  geom_point(aes(fill = name), size = 3, shape = 21) + #
  # setting colours
  scale_color_manual(values = c("Microbe" = "#009E73",  "Plant" = "#F0E422",  "Invert"= "#0072B2",  "Vert" = "#D55E00" )) +
  scale_fill_manual(values = c("Microbe" = "#009E73",  "Plant" = "#F0E422",  "Invert"= "#0072B2",  "Vert" = "#D55E00")) +
  scale_y_discrete(labels = c("Microbe" = "Microbe",  "Plant" = "Plant",  "Invert"= "Invertebrate",  "Vert" = "Vertebrate")) +
  annotate('text', x = 0.93, y = 1:4 + 0.15, label= paste("italic(k)==", res_symbiont_tax_broad1$n), parse=TRUE, hjust = "left", size=3.5) +
  labs(x = expression(paste(italic(r), " (correlation)")), y = "", size = expression(paste(italic(n), " (# of species)")), tag = "d") +
  guides(fill = "none", colour = "none") +
  theme_bw() +
  theme(legend.position="none") +
  theme(axis.text.y = element_text(size = 10, colour ="black",hjust = 0.5, angle = 90)) +
  # putting pictures in
  annotation_custom(rasterGrob(image_microbe_parasite), xmin = -1, xmax = -0.8, ymin = 0.6, ymax = 1.2) + 
  annotation_custom(rasterGrob(image_plant_parasite), xmin = -1, xmax = -0.8, ymin = 1.6, ymax = 2.2) +
  annotation_custom(rasterGrob(image_invertebrate_parasite), xmin = -1, xmax = -0.8, ymin = 2.6, ymax = 3.2) + 
  annotation_custom(rasterGrob(image_vertebrate_parasite), xmin = -1, xmax = -0.8, ymin = 3.6, ymax = 4.2)

Figure 2d: A forest plot showing the group-wise means (the categorical variable symbiont_tax_broad) with their 95% confidences interval (thick lines) and 95% prediction intervals (thin lines) with observed effect sizes based on various sample sizes.

Testing specialization 1: host range

# meta-regression
mr_host_range_link_ratio <- rma.mv(yi = Zr, V = VZr, mods = ~log(host_range_link_ratio), 
    random = ~1 | authors, data = dat)

Supplementary Table 5: Regression coefficients (Estimate), 95% confidence intervals (CIs), variance components (V) and variance explained, R²_[marginal] (R2) from the meta-regression with log(host_range_link_ratio).

# getting marginal R2
r2_host_range_link_ratio <- R2(mr_host_range_link_ratio)

# getting estimates: name does not work for slopes
res_host_range_link_ratio <- get_est(mr_host_range_link_ratio, mod = "log(host_range_link_ratio)")

# creating a table
tibble(`Fixed effect` = c("Intercept", "log(host_range_link_ratio)"), Estimate = c(res_host_range_link_ratio$estimate), 
    `Lower CI [0.025]` = c(res_host_range_link_ratio$lowerCL), `Upper CI  [0.975]` = c(res_host_range_link_ratio$upperCL), 
    `V[authors]` = c(mr_host_range_link_ratio$sigma2, NA), R2 = c(r2_host_range_link_ratio[1], 
        NA)) %>% kable("html", digits = 3) %>% kable_styling("striped", position = "left")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	V[authors]	R2
Intercept	0.360	0.318	0.401	0.033	0.012
log(host_range_link_ratio)	-0.057	-0.150	0.037	NA	NA

# newmods <- seq(-0.3, 2.2, by = 0.1)
# pred_host_range_link_ratio <-predict.rma(mr_host_range_link_ratio, newmods = newmods) 
# ribbon_dat <- tibble(newmods = newmods, ymin = pred_host_range_link_ratio$ci.lb, ymax = pred_host_range_link_ratio$ci.ub)
pred_host_range_link_ratio <-predict.rma(mr_host_range_link_ratio) 

# plotting

fig_host_range_link_ratio <-  dat %>% 
  filter(!is.na(host_range_link_ratio))  %>% # getting ride of NA values
  mutate(ymin = pred_host_range_link_ratio$ci.lb, 
         ymax = pred_host_range_link_ratio$ci.ub,
         ymin2 = pred_host_range_link_ratio$cr.lb,
         ymax2 = pred_host_range_link_ratio$cr.ub,
         pred = pred_host_range_link_ratio$pred) %>% 
  ggplot(aes(x = log(host_range_link_ratio), y = Zr, size = (1/VZr) + 3, )) +
  geom_point(shape = 21, fill = "grey90") +
  #geom_ribbon(aes(ymin = ymin, ymax = ymax), fill = "#0072B2")  + # not quite sure why this does not work
  geom_smooth(aes(y = ymin2), method =  "loess", se = FALSE, lty =  "dotted", lwd = 0.25, colour = "#0072B2") +
  geom_smooth(aes(y = ymax2), method =  "loess", se = FALSE, lty = "dotted", lwd = 0.25, colour = "#0072B2") +
  geom_smooth(aes(y = ymin), method =  "loess", se = FALSE,lty = "dotted", lwd = 0.25, colour ="#D55E00") +
  geom_smooth(aes(y = ymax), method =  "loess", se = FALSE, lty ="dotted", lwd = 0.25, colour ="#D55E00") + 
  geom_smooth(aes(y = pred), method =  "loess", se = FALSE, lty ="dashed", lwd = 0.5, colour ="black") +  
  ylim(-1, 2) + xlim(-0.05, 2) +
  #geom_abline(intercept = mr_host_range_link_ratio$beta[[1]], slope = mr_host_range_link_ratio$beta[[2]], alpha = 0.7, linetype = "dashed", size = 0.5) +
  labs(x = "ln(host range link ratio)", y = expression(paste(italic(Zr), " (effect size)")), size = expression(paste(italic(n), " (# of species)"))) +
  guides(fill = "none", colour = "none") +
  # themses
  theme_bw() +
  theme(legend.position= c(1, 1), legend.justification = c(1, 1)) +
  theme(legend.direction="horizontal") +
  #theme(legend.background = element_rect(fill = "white", colour = "black")) +
  theme(legend.background = element_blank()) +
  theme(axis.text.y = element_text(size = 10, colour ="black", hjust = 0.5, angle = 90)) 

fig_host_range_link_ratio

Supplementary Figure 1: A bubble plot showing a predicted regression line for the contentious variable log(host_range_link_ratio) with their 95% confidences regions (orange dotted lines) and 95% prediction regions (blue dotted lines) with observed effect sizes based on various sample sizes.

Testing specialization 2: taxonomic breadth

# meta-regression
mr_host_range_taxonomic_breadth <- rma.mv(yi = Zr, V = VZr, mods = ~log(host_range_taxonomic_breadth), 
    random = ~1 | authors, data = dat)

Supplementary Table 6: Regression coefficients (Estimate), 95% confidence intervals (CIs), variance components (V) and variance explained, R²_[marginal] (R2) from the meta-regression with log(host_range_taxonomic_breadth).

# getting marginal R2
r2_host_range_taxonomic_breadth <- R2(mr_host_range_taxonomic_breadth)

# getting estimates: name does not work for slopes
res_host_range_taxonomic_breadth <- get_est(mr_host_range_taxonomic_breadth, 
    mod = "log(host_range_taxonomic_breadth)")

# creating a table
tibble(`Fixed effect` = c("Intercept", "log(host_range_taxonomic_breadth)"), 
    Estimate = c(res_host_range_taxonomic_breadth$estimate), `Lower CI [0.025]` = c(res_host_range_taxonomic_breadth$lowerCL), 
    `Upper CI  [0.975]` = c(res_host_range_taxonomic_breadth$upperCL), `V[authors]` = c(mr_host_range_taxonomic_breadth$sigma2, 
        NA), R2 = c(r2_host_range_taxonomic_breadth[1], NA)) %>% kable("html", 
    digits = 3) %>% kable_styling("striped", position = "left")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	V[authors]	R2
Intercept	0.364	0.323	0.405	0.032	0.008
log(host_range_taxonomic_breadth)	-0.056	-0.170	0.057	NA	NA

pred_host_range_taxonomic_breadth <-predict.rma(mr_host_range_taxonomic_breadth) 

# plotting

fig_host_range_taxonomic_breadth <-  dat %>% 
  filter(!is.na(host_range_taxonomic_breadth))  %>% # getting ride of NA values
  mutate(ymin = pred_host_range_taxonomic_breadth$ci.lb, 
         ymax = pred_host_range_taxonomic_breadth$ci.ub,
         ymin2 = pred_host_range_taxonomic_breadth$cr.lb,
         ymax2 = pred_host_range_taxonomic_breadth$cr.ub,
         pred = pred_host_range_taxonomic_breadth$pred) %>% 
  ggplot(aes(x = log(host_range_taxonomic_breadth), y = Zr, size = (1/VZr) + 3, )) +
  geom_point(shape = 21, fill = "grey90") +
  #geom_ribbon(aes(ymin = ymin, ymax = ymax), fill = "#0072B2")  + # not quite sure why this does not work
  geom_smooth(aes(y = ymin2), method =  "loess", se = FALSE, lty =  "dotted", lwd = 0.25, colour = "#0072B2") +
  geom_smooth(aes(y = ymax2), method =  "loess", se = FALSE, lty = "dotted", lwd = 0.25, colour = "#0072B2") +
  geom_smooth(aes(y = ymin), method =  "loess", se = FALSE,lty = "dotted", lwd = 0.25, colour ="#D55E00") +
  geom_smooth(aes(y = ymax), method =  "loess", se = FALSE, lty ="dotted", lwd = 0.25, colour ="#D55E00") + 
  geom_smooth(aes(y = pred), method =  "loess", se = FALSE, lty ="dashed", lwd = 0.5, colour ="black") +  
  ylim(-1, 2) + xlim(0, 1.5) +
  #geom_abline(intercept = mr_host_range_link_ratio$beta[[1]], slope = mr_host_range_link_ratio$beta[[2]], alpha = 0.7, linetype = "dashed", size = 0.5) +
  labs(x = "ln(host range taxonomic breadth)", y = expression(paste(italic(Zr), " (effect size)")), size = expression(paste(italic(n), " (# of species)"))) +
  guides(fill = "none", colour = "none") +
  # themses
  theme_bw() +
  theme(legend.position= c(1, 1), legend.justification = c(1, 1)) +
  theme(legend.direction="horizontal") +
  #theme(legend.background = element_rect(fill = "white", colour = "black")) +
  theme(legend.background = element_blank()) +
  theme(axis.text.y = element_text(size = 10, colour ="black", hjust = 0.5, angle = 90)) 

fig_host_range_taxonomic_breadth

Supplementary Figure 2: A bubble plot showing a predicted regression line for the contentious variable log(log(host_range_taxonomic_breadth) with their 95% confidences regions (orange dotted lines) and 95% prediction regions (blue dotted lines) with observed effect sizes based on various sample sizes.

The place of symbiosis: ednosymbiosis vs. ectosymbiosis

# reordering
dat$endo_or_ecto <- factor(dat$endo_or_ecto, levels = c("Endo/Ecto", "Endo", 
    "Ecto"))

# meta-regression: mutiple intercepts
mr_endo_or_ecto1 <- rma.mv(yi = Zr, V = VZr, mods = ~endo_or_ecto - 1, test = "t", 
    random = ~1 | authors, data = dat)

# meta-regression: contrast 1
mr_endo_or_ecto2 <- rma.mv(yi = Zr, V = VZr, mods = ~endo_or_ecto, test = "t", 
    random = ~1 | authors, data = dat)

# meta-regression: contrast 2
mr_endo_or_ecto3 <- rma.mv(yi = Zr, V = VZr, mods = ~relevel(endo_or_ecto, ref = "Endo"), 
    test = "t", random = ~1 | authors, data = dat)

Supplementary Table 7: Regression coefficients (Estimate), 95% confidence intervals (CIs), variance components (V) and variance explained, R²_[marginal] (R2) from the meta-regression with endo_or_ecto.

# getting marginal R2
r2_endo_or_ecto1 <- R2(mr_endo_or_ecto1)

# getting estimates
res_endo_or_ecto1 <- get_est(mr_endo_or_ecto1, mod = "endo_or_ecto")
res_endo_or_ecto2 <- get_est(mr_endo_or_ecto2, mod = "endo_or_ecto")
res_endo_or_ecto3 <- get_est(mr_endo_or_ecto3, mod = "endo_or_ecto")

# creating a table
tibble(`Fixed effect` = c(as.character(res_endo_or_ecto1$name), cont_gen(res_endo_or_ecto1$name)), 
    Estimate = c(res_endo_or_ecto1$estimate, res_endo_or_ecto2$estimate[-1], 
        res_endo_or_ecto3$estimate[-(1:2)]), `Lower CI [0.025]` = c(res_endo_or_ecto1$lowerCL, 
        res_endo_or_ecto2$lowerCL[-1], res_endo_or_ecto3$lowerCL[-(1:2)]), `Upper CI  [0.975]` = c(res_endo_or_ecto1$upperCL, 
        res_endo_or_ecto2$upperCL[-1], res_endo_or_ecto3$upperCL[-(1:2)]), `V[authors]` = c(mr_endo_or_ecto1$sigma2, 
        rep(NA, 5)), R2 = c(r2_endo_or_ecto1[1], rep(NA, 5))) %>% kable("html", 
    digits = 3) %>% kable_styling("striped", position = "left")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	V[authors]	R2
Endo/Ecto	0.334	0.035	0.632	0.034	0.001
Endo	0.342	0.298	0.387	NA	NA
Ecto	0.353	0.287	0.418	NA	NA
Endo/Ecto-Endo	0.009	-0.293	0.311	NA	NA
Endo/Ecto-Ecto	0.019	-0.287	0.325	NA	NA
Endo-Ecto	0.010	-0.069	0.090	NA	NA

# getting images
image_endoparasite <- readPNG("../images/endoparasite_transparentbg.png")
image_ectoparasite <- readPNG("../images/ectoparasite_transparentbg.png")

# adding sample size (k) for each category
k_endo_or_ecto <- dat %>% group_by(endo_or_ecto) %>% count()
# getting estimates and predicitons
pred_endo_or_ecto <- get_pred(mr_endo_or_ecto1, mod = "endo_or_ecto") 
res_endo_or_ecto1 <- left_join(res_endo_or_ecto1, k_endo_or_ecto, by =  c("name" = "endo_or_ecto"))  %>% left_join(pred_endo_or_ecto)
#res_symbiosis1 
# drawing a funnel plot - fig 2b
fig_endo_or_ecto <- ggplot(data = res_endo_or_ecto1, aes(x = tanh(estimate), y = name)) +
  scale_x_continuous(limits=c(-1, 1), breaks = seq(-1, 1, by = 0.2) ) +
  geom_quasirandom(data = dat %>% filter(!is.na(endo_or_ecto)), 
                   aes(x= tanh(Zr), y = endo_or_ecto, size = ((1/VZr) + 3), colour = endo_or_ecto), groupOnX = FALSE, alpha=0.4) + 
  # 95 %precition interval (PI)
  geom_errorbarh(aes(xmin = tanh(lowerPR), xmax = tanh(upperPR)),  height = 0, show.legend = F, size = 0.5, alpha = 0.6) +
  # 95 %CI
  geom_errorbarh(aes(xmin = tanh(lowerCL), xmax = tanh(upperCL)),  height = 0, show.legend = F, size = 1.2) +
  geom_vline(xintercept = 0, linetype = 2, colour = "black", alpha = 0.3) +
  # creating dots and different size (bee-swarm and bubbles)
  geom_point(aes(fill = name), size = 3, shape = 21) + #
  # setting colours
  scale_color_manual(values = c("Endo/Ecto" = "#0072B2",  "Endo" = "#D55E00",  "Ecto"= "#CC79A7")) +
  scale_fill_manual(values = c("Endo/Ecto" = "#0072B2",  "Endo" = "#D55E00",  "Ecto"= "#CC79A7")) +
  scale_y_discrete(labels = c("Endo/Ecto" = "Both",  "Endo" = "Endosymbiosis",  "Ecto"= "Ectosymbiosis")) +
  annotate('text', x = 0.93, y = 1:3 + 0.15, label= paste("italic(k)==", res_endo_or_ecto1$n), parse=TRUE, hjust = "left", size=3.5) +
  labs(x = expression(paste(italic(r), " (correlation)")), y = "", size = expression(paste(italic(n), " (# of species)")) ) +
  guides(fill = "none", colour = "none") +
  theme_bw() +
  theme(legend.position= c(0, 1), legend.justification = c(0,1)) +
  theme(legend.direction = "horizontal") +
  #theme(legend.background = element_rect(fill = "white", colour = "black")) +
  theme(legend.background = element_blank()) +
  theme(axis.text.y = element_text(size = 10, colour ="black", hjust = 0.5, angle = 90)) +
  # adding images
  annotation_custom(rasterGrob(image_endoparasite), xmin = -1, xmax = -0.8, ymin = 1.6, ymax = 2.2) +
  annotation_custom(rasterGrob(image_ectoparasite), xmin = -1, xmax = -0.8, ymin = 2.6, ymax = 3.2) + 
  annotation_custom(rasterGrob(image_ectoparasite), xmin = -1.1, xmax = -0.9, ymin = 0.6, ymax = 1.2) + 
  annotation_custom(rasterGrob(image_endoparasite), xmin = -0.9, xmax = -0.7, ymin = 0.6, ymax = 1.2)

fig_endo_or_ecto

# fig 3e

e <- ggplot(data = res_endo_or_ecto1, aes(x = tanh(estimate), y = name)) +
  scale_x_continuous(limits=c(-1, 1), breaks = seq(-1, 1, by = 0.2) ) +
  geom_quasirandom(data = dat %>% filter(!is.na(endo_or_ecto)), 
                   aes(x= tanh(Zr), y = endo_or_ecto, size = ((1/VZr) + 3), colour = endo_or_ecto), groupOnX = FALSE, alpha=0.4) + 
  # 95 %precition interval (PI)
  geom_errorbarh(aes(xmin = tanh(lowerPR), xmax = tanh(upperPR)),  height = 0, show.legend = F, size = 0.5, alpha = 0.6) +
  # 95 %CI
  geom_errorbarh(aes(xmin = tanh(lowerCL), xmax = tanh(upperCL)),  height = 0, show.legend = F, size = 1.2) +
  geom_vline(xintercept = 0, linetype = 2, colour = "black", alpha = 0.3) +
  # creating dots and different size (bee-swarm and bubbles)
  geom_point(aes(fill = name), size = 3, shape = 21) + #
  # setting colours
  scale_color_manual(values = c("Endo/Ecto" = "#0072B2",  "Endo" = "#D55E00",  "Ecto"= "#CC79A7")) +
  scale_fill_manual(values = c("Endo/Ecto" = "#0072B2",  "Endo" = "#D55E00",  "Ecto"= "#CC79A7")) +
  scale_y_discrete(labels = c("Endo/Ecto" = "Both",  "Endo" = "Endosymbiosis",  "Ecto"= "Ectosymbiosis")) +
  annotate('text', x = 0.93, y = 1:3 + 0.15, label= paste("italic(k)==", res_endo_or_ecto1$n), parse=TRUE, hjust = "left", size=3.5) +
  labs(x = "", y = "", size = expression(paste(italic(n), " (# of species)")), tag = "e" ) +
  guides(fill = "none", colour = "none") +
  theme_bw() +
  theme(legend.position="none") +
  theme(axis.text.y = element_text(size = 10, colour ="black",hjust = 0.5, angle = 90)) +
  # adding images
  annotation_custom(rasterGrob(image_endoparasite), xmin = -1, xmax = -0.8, ymin = 1.6, ymax = 2.2) +
  annotation_custom(rasterGrob(image_ectoparasite), xmin = -1, xmax = -0.8, ymin = 2.6, ymax = 3.2) + 
  annotation_custom(rasterGrob(image_ectoparasite), xmin = -1.1, xmax = -0.9, ymin = 0.6, ymax = 1.2) + 
  annotation_custom(rasterGrob(image_endoparasite), xmin = -0.9, xmax = -0.7, ymin = 0.6, ymax = 1.2)

Figure 3e: A forest plot showing the group-wise means (the categorical variable endo_or_ecto) with their 95% confidences interval (thick lines) and 95% prediction intervals (thin lines) with observed effect sizes based on various sample sizes.

The effect of the model of transmission

# meta-regression: mutiple intercepts
mr_mode_of_transmission_broad1 <- rma.mv(yi = Zr, V = VZr, mods = ~mode_of_transmission_broad - 
    1, test = "t", random = ~1 | authors, data = dat)

# meta-regression: contrast 1
mr_mode_of_transmission_broad2 <- rma.mv(yi = Zr, V = VZr, mods = ~mode_of_transmission_broad, 
    test = "t", random = ~1 | authors, data = dat)

# meta-regression: contrast 2
mr_mode_of_transmission_broad3 <- rma.mv(yi = Zr, V = VZr, mods = ~relevel(mode_of_transmission_broad, 
    ref = "vertical"), test = "t", random = ~1 | authors, data = dat)

Supplementary Table 8: Regression coefficients (Estimate), 95% confidence intervals (CIs), variance components (V) and variance explained, R²_[marginal] (R2) from the meta-regression with mode_of_transmission_broad.

# getting marginal R2
r2_mode_of_transmission_broad1 <- R2(mr_mode_of_transmission_broad1)

# getting estimates
res_mode_of_transmission_broad1 <- get_est(mr_mode_of_transmission_broad1, mod = "mode_of_transmission_broad")
res_mode_of_transmission_broad2 <- get_est(mr_mode_of_transmission_broad2, mod = "mode_of_transmission_broad")
res_mode_of_transmission_broad3 <- get_est(mr_mode_of_transmission_broad3, mod = "mode_of_transmission_broad")

# creating a table
tibble(`Fixed effect` = c(as.character(res_mode_of_transmission_broad1$name), 
    cont_gen(res_mode_of_transmission_broad1$name)), Estimate = c(res_mode_of_transmission_broad1$estimate, 
    res_mode_of_transmission_broad2$estimate[-1], res_mode_of_transmission_broad3$estimate[-(1:2)]), 
    `Lower CI [0.025]` = c(res_mode_of_transmission_broad1$lowerCL, res_mode_of_transmission_broad2$lowerCL[-1], 
        res_mode_of_transmission_broad3$lowerCL[-(1:2)]), `Upper CI  [0.975]` = c(res_mode_of_transmission_broad1$upperCL, 
        res_mode_of_transmission_broad2$upperCL[-1], res_mode_of_transmission_broad3$upperCL[-(1:2)]), 
    `V[authors]` = c(mr_mode_of_transmission_broad1$sigma2, rep(NA, 5)), R2 = c(r2_mode_of_transmission_broad1[1], 
        rep(NA, 5))) %>% kable("html", digits = 3) %>% kable_styling("striped", 
    position = "left")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	V[authors]	R2
both	0.371	0.298	0.443	0.024	0.219
horizontal	0.285	0.242	0.328	NA	NA
vertical	0.496	0.421	0.570	NA	NA
both-horizontal	-0.085	-0.170	-0.001	NA	NA
both-vertical	0.125	0.021	0.229	NA	NA
horizontal-vertical	-0.210	-0.297	-0.124	NA	NA

# getting images
image_horizontal <- readPNG("../images/horizontal_transparentbg.png")
image_vertical <- readPNG("../images/vertical_transparentbg.png")
image_both <- readPNG("../images/horizontal_vertical_transparentbg.png")
# adding sample size (k) for each category
k_mode_of_transmission_broad <- dat %>% group_by(mode_of_transmission_broad) %>% count()
# getting estimates and predicitons
pred_mode_of_transmission_broad <- get_pred(mr_mode_of_transmission_broad1, mod = "mode_of_transmission_broad") 
res_mode_of_transmission_broad1 <- left_join(res_mode_of_transmission_broad1, k_mode_of_transmission_broad, by =  c("name" = "mode_of_transmission_broad"))  %>% left_join(pred_mode_of_transmission_broad)
#res_symbiosis1 
# drawing a funnel plot - fig 2b
fig_mode_of_transmission_broad <- ggplot(data = res_mode_of_transmission_broad1, aes(x = tanh(estimate), y = name)) +
  scale_x_continuous(limits=c(-1, 1), breaks = seq(-1, 1, by = 0.2) ) +
  geom_quasirandom(data = dat %>% filter(!is.na(mode_of_transmission_broad)), 
                   aes(x= tanh(Zr), y = mode_of_transmission_broad, size = ((1/VZr) + 3), colour = mode_of_transmission_broad), groupOnX = FALSE, alpha=0.4) + 
  # 95 %precition interval (PI)
  geom_errorbarh(aes(xmin = tanh(lowerPR), xmax = tanh(upperPR)),  height = 0, show.legend = F, size = 0.5, alpha = 0.6) +
  # 95 %CI
  geom_errorbarh(aes(xmin = tanh(lowerCL), xmax = tanh(upperCL)),  height = 0, show.legend = F, size = 1.2) +
  geom_vline(xintercept = 0, linetype = 2, colour = "black", alpha = 0.3) +
  # creating dots and different size (bee-swarm and bubbles)
  geom_point(aes(fill = name), size = 3, shape = 21) + #
  # setting colours
  scale_color_manual(values = c("both" = "#0072B2",  "horizontal" = "#D55E00",  "vertical"= "#CC79A7")) +
  scale_fill_manual(values = c("both" = "#0072B2",  "horizontal" = "#D55E00",  "vertical"= "#CC79A7")) +
  scale_y_discrete(labels = c("both" = "Both",  "horizontal" = "Horizontal",  "vertical"= "Vertical")) +
  annotate('text', x = 0.93, y = (1:3 + 0.15), label= paste("italic(k)==", res_mode_of_transmission_broad1$n), parse=TRUE, hjust = "left", size=3.5) +
  labs(x = expression(paste(italic(r), " (correlation)")), y = "", size = expression(paste(italic(n), " (# of species)")) ) +
  guides(fill = "none", colour = "none") +
  theme_bw() +
  theme(legend.position= c(0, 1), legend.justification = c(0,1)) +
  theme(legend.direction = "horizontal") +
  #theme(legend.background = element_rect(fill = "white", colour = "black")) +
  theme(legend.background = element_blank()) +
  theme(axis.text.y = element_text(size = 10, colour ="black", hjust = 0.5, angle = 90)) +
  # adding images
  annotation_custom(rasterGrob(image_horizontal), xmin = -1, xmax = -0.7, ymin = 1.4, ymax = 2.2) +
  annotation_custom(rasterGrob(image_vertical), xmin = -1, xmax = -0.7, ymin = 2.4, ymax = 3.2) + 
  annotation_custom(rasterGrob(image_both), xmin = -1, xmax = -0.7, ymin = 0.4, ymax = 1.2) 

fig_mode_of_transmission_broad

# fig 3f

f <- ggplot(data = res_mode_of_transmission_broad1, aes(x = tanh(estimate), y = name)) +
  scale_x_continuous(limits=c(-1, 1), breaks = seq(-1, 1, by = 0.2) ) +
  geom_quasirandom(data = dat %>% filter(!is.na(mode_of_transmission_broad)), 
                   aes(x= tanh(Zr), y = mode_of_transmission_broad, size = ((1/VZr) + 3), colour = mode_of_transmission_broad), groupOnX = FALSE, alpha=0.4) + 
  # 95 %precition interval (PI)
  geom_errorbarh(aes(xmin = tanh(lowerPR), xmax = tanh(upperPR)),  height = 0, show.legend = F, size = 0.5, alpha = 0.6) +
  # 95 %CI
  geom_errorbarh(aes(xmin = tanh(lowerCL), xmax = tanh(upperCL)),  height = 0, show.legend = F, size = 1.2) +
  geom_vline(xintercept = 0, linetype = 2, colour = "black", alpha = 0.3) +
  # creating dots and different size (bee-swarm and bubbles)
  geom_point(aes(fill = name), size = 3, shape = 21) + #
  # setting colours
  scale_color_manual(values = c("both" = "#0072B2",  "horizontal" = "#D55E00",  "vertical"= "#CC79A7")) +
  scale_fill_manual(values = c("both" = "#0072B2",  "horizontal" = "#D55E00",  "vertical"= "#CC79A7")) +
  scale_y_discrete(labels = c("both" = "Both",  "horizontal" = "Horizontal",  "vertical"= "Vertical")) +
  annotate('text', x = 0.93, y = (1:3 + 0.15), label= paste("italic(k)==", res_mode_of_transmission_broad1$n), parse=TRUE, hjust = "left", size=3.5) +
  labs(x = "", y = "", size = expression(paste(italic(n), " (# of species)")), tag = "f" ) +
  guides(fill = "none", colour = "none") +
  theme_bw() +
  theme(legend.position="none") +
  theme(axis.text.y = element_text(size = 10, colour ="black",hjust = 0.5, angle = 90)) +
  # adding images
  annotation_custom(rasterGrob(image_horizontal), xmin = -1, xmax = -0.7, ymin = 1.4, ymax = 2.2) +
  annotation_custom(rasterGrob(image_vertical), xmin = -1, xmax = -0.7, ymin = 2.4, ymax = 3.2) + 
  annotation_custom(rasterGrob(image_both), xmin = -1, xmax = -0.7, ymin = 0.4, ymax = 1.2) 

Figure 2f: A forest plot showing the group-wise means (the categorical variable mode_of_transmission_broad) with their 95% confidences interval (thick lines) and 95% prediction intervals (thin lines) with observed effect sizes based on various sample sizes.

The combined effect of symbiosis and mode of transmission

# reordering
dat$symbiosis_transmission <- factor(dat$symbiosis_transmission, levels = c("Mutualistboth", 
    "Mutualisthorizontal", "Mutualistvertical", "Parasiteboth", "Parasitehorizontal"), 
    labels = c("MutualistBoth", "MutualistHorizontal", "MutualistVertical", 
        "ParasiteBoth", "ParasiteHorizontal"))

# meta-regression: mutiple intercepts
mr_symbiosis_transmission1 <- rma.mv(yi = Zr, V = VZr, mods = ~symbiosis_transmission - 
    1, test = "t", random = ~1 | authors, data = dat)

# # meta-regression: contrasts x 10 getting the level names out
level_names <- levels(dat$symbiosis_transmission)

# helper function to run metafor meta-regression
run_rma <- function(name) {
    rma.mv(yi = Zr, V = VZr, mods = ~relevel(symbiosis_transmission, ref = name), 
        test = "t", random = ~1 | authors, data = dat)
}

# results of meta-regression including all contrast results; taking the last
# level out ([-length(level_names)])
mr_symbiosis_transmission <- map(level_names[-length(level_names)], run_rma)

Supplementary Table 9: Regression coefficients (Estimate), 95% confidence intervals (CIs), variance components (V) and variance explained, R²_[marginal] (R2) from the meta-regression with symbiosis_transmission.

# getting marginal R2
r2_symbiosis_transmission1 <- R2(mr_symbiosis_transmission1)

# getting estimates
res_symbiosis_transmission1 <- get_est(mr_symbiosis_transmission1, mod = "symbiosis_transmission")
res_symbiosis_transmission <- map(mr_symbiosis_transmission, ~get_est(.x, mod = "symbiosis_transmission"))

# a list of the numbers to take out unnecessary contrasts
contra_list <- Map(seq, from = 1, to = 1:4)

# you need to flatten twice: first to make it a list and make it a vector
estimates <- map2(res_symbiosis_transmission, contra_list, ~.x[-(.y), "estimate"]) %>% 
    flatten() %>% flatten_dbl()
lowerCLs <- map2(res_symbiosis_transmission, contra_list, ~.x[-(.y), "lowerCL"]) %>% 
    flatten() %>% flatten_dbl()
upperCLs <- map2(res_symbiosis_transmission, contra_list, ~.x[-(.y), "upperCL"]) %>% 
    flatten() %>% flatten_dbl()

# creating a table
tibble(`Fixed effect` = c(as.character(res_symbiosis_transmission1$name), cont_gen(res_symbiosis_transmission1$name)), 
    Estimate = c(res_symbiosis_transmission1$estimate, estimates), `Lower CI [0.025]` = c(res_symbiosis_transmission1$lowerCL, 
        lowerCLs), `Upper CI  [0.975]` = c(res_symbiosis_transmission1$upperCL, 
        upperCLs), `V[authors]` = c(mr_symbiosis_transmission1$sigma2, rep(NA, 
        (5 + choose(5, 2)) - 1)), R2 = c(r2_symbiosis_transmission1[1], rep(NA, 
        (5 + choose(5, 2)) - 1))) %>% kable("html", digits = 3) %>% kable_styling("striped", 
    position = "left") %>% scroll_box(width = "100%", height = "300px")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	V[authors]	R2
MutualistBoth	0.453	0.224	0.682	0.024	0.213
MutualistHorizontal	0.290	0.189	0.391	NA	NA
MutualistVertical	0.493	0.417	0.568	NA	NA
ParasiteBoth	0.362	0.285	0.439	NA	NA
ParasiteHorizontal	0.284	0.238	0.331	NA	NA
MutualistBoth-MutualistHorizontal	-0.163	-0.413	0.088	NA	NA
MutualistBoth-MutualistVertical	0.040	-0.201	0.281	NA	NA
MutualistBoth-ParasiteBoth	-0.091	-0.333	0.150	NA	NA
MutualistBoth-ParasiteHorizontal	-0.168	-0.402	0.065	NA	NA
MutualistHorizontal-MutualistVertical	0.202	0.076	0.329	NA	NA
MutualistHorizontal-ParasiteBoth	0.071	-0.056	0.198	NA	NA
MutualistHorizontal-ParasiteHorizontal	-0.006	-0.115	0.103	NA	NA
MutualistVertical-ParasiteBoth	-0.131	-0.239	-0.023	NA	NA
MutualistVertical-ParasiteHorizontal	-0.208	-0.297	-0.119	NA	NA
ParasiteBoth-ParasiteHorizontal	-0.077	-0.167	0.013	NA	NA

# colour list
colour_ls <- c("#000000", "#E69F00", "#56B4E9", "#009E73",  "#F0E422",  "#0072B2",  "#D55E00", "#CC79A7", "#00008B", "#8B0A50", "#54FF9F", "#999999")

# adding sample size (k) for each category
k_symbiosis_transmission <- dat %>% group_by(symbiosis_transmission) %>% count()
# getting estimates and predicitons
pred_symbiosis_transmission <- get_pred(mr_symbiosis_transmission1, mod = "symbiosis_transmission") 
res_symbiosis_transmission1 <- left_join(res_symbiosis_transmission1, k_symbiosis_transmission, by =  c("name" = "symbiosis_transmission"))  %>% left_join(pred_symbiosis_transmission)
#res_symbiosis1 
# drawing a funnel plot - fig 2b
fig_symbiosis_transmission <- ggplot(data = res_symbiosis_transmission1, aes(x = tanh(estimate), y = name)) +
  scale_x_continuous(limits=c(-1, 1), breaks = seq(-1, 1, by = 0.2) ) +
  geom_quasirandom(data = dat %>% filter(!is.na(symbiosis_transmission)), 
                   aes(x= tanh(Zr), y = symbiosis_transmission, size = ((1/VZr) + 3), colour = symbiosis_transmission), groupOnX = FALSE, alpha=0.4) + 
  # 95 %precition interval (PI)
  geom_errorbarh(aes(xmin = tanh(lowerPR), xmax = tanh(upperPR)),  height = 0, show.legend = F, size = 0.5, alpha = 0.6) +
  # 95 %CI
  geom_errorbarh(aes(xmin = tanh(lowerCL), xmax = tanh(upperCL)),  height = 0, show.legend = F, size = 1.2) +
  geom_vline(xintercept = 0, linetype = 2, colour = "black", alpha = 0.3) +
  # creating dots and different size (bee-swarm and bubbles)
  geom_point(aes(fill = name), size = 3, shape = 21) + #
  # setting colours
  scale_color_manual(values =   c("MutualistBoth"= colour_ls[1], "MutualistHorizontal"= colour_ls[2], "MutualistVertical" = colour_ls[3],"ParasiteBoth"= colour_ls[4], "ParasiteHorizontal" = colour_ls[5])) +
  scale_fill_manual(values = c("MutualistBoth"= colour_ls[1], "MutualistHorizontal"= colour_ls[2], "MutualistVertical" = colour_ls[3],"ParasiteBoth"= colour_ls[4], "ParasiteHorizontal" = colour_ls[5])) +
  scale_y_discrete(labels = c("MutualistBoth" = "Mutualist-\nBoth", "MutualistHorizontal" = "Mutualist-\nHorizontal","MutualistVertical" = "Mutualist-\nVertical", "ParasiteBoth" = "Parasite-\nBoth", "ParasiteHorizontal" = "Parasite-\nHorizontal")) +
  annotate('text', x = 0.93, y = 1:5 + 0.15, label= paste("italic(k)==", res_symbiosis_transmission1$n), parse=TRUE, hjust = "left", size=3.5) +
  labs(x = expression(paste(italic(r), " (correlation)")), y = "", size = expression(paste(italic(n), " (# of species)")) ) +
  guides(fill = "none", colour = "none") +
  theme_bw() +
  theme(legend.position= c(0, 1), legend.justification = c(0,1)) +
  theme(legend.direction="horizontal") +
  #theme(legend.background = element_rect(fill = "white", colour = "black")) +
  theme(legend.background = element_blank()) +
  theme(axis.text.y = element_text(size = 10, colour ="black", hjust = 0.5, angle = 90)) +
  # putting pictures in
  annotation_custom(rasterGrob(image_mutualism), xmin = -1.1, xmax = -0.9, ymin = 0.6, ymax = 1.2) + 
  annotation_custom(rasterGrob(image_both), xmin = -0.9, xmax = -0.6, ymin = 0.4, ymax = 1.2) + 
  annotation_custom(rasterGrob(image_mutualism), xmin = -1.1, xmax = -0.9, ymin = 1.6, ymax = 2.2) +
  annotation_custom(rasterGrob(image_horizontal),xmin = -0.9, xmax = -0.6, ymin = 1.4, ymax = 2.2) +
  annotation_custom(rasterGrob(image_mutualism), xmin = -1.1, xmax = -0.9, ymin = 2.6, ymax = 3.2) + 
  annotation_custom(rasterGrob(image_vertical), xmin = -0.9, xmax = -0.6, ymin = 2.4, ymax = 3.2) + 
  annotation_custom(rasterGrob(image_parasitism), xmin = -1.1, xmax = -0.9, ymin = 3.6, ymax = 4.2) +
  annotation_custom(rasterGrob(image_both), xmin = -0.9, xmax = -0.6, ymin = 3.4, ymax = 4.2) +
  annotation_custom(rasterGrob(image_parasitism), xmin = -1.1, xmax = -0.9, ymin = 4.6, ymax = 5.2) +
  annotation_custom(rasterGrob(image_horizontal), xmin = -0.9, xmax = -0.6, ymin = 4.4, ymax = 5.2)

fig_symbiosis_transmission

## fig 3
g <- ggplot(data = res_symbiosis_transmission1, aes(x = tanh(estimate), y = name)) +
  scale_x_continuous(limits=c(-1, 1), breaks = seq(-1, 1, by = 0.2) ) +
  geom_quasirandom(data = dat %>% filter(!is.na(symbiosis_transmission)), 
                   aes(x= tanh(Zr), y = symbiosis_transmission, size = ((1/VZr) + 3), colour = symbiosis_transmission), groupOnX = FALSE, alpha=0.4) + 
  # 95 %precition interval (PI)
  geom_errorbarh(aes(xmin = tanh(lowerPR), xmax = tanh(upperPR)),  height = 0, show.legend = F, size = 0.5, alpha = 0.6) +
  # 95 %CI
  geom_errorbarh(aes(xmin = tanh(lowerCL), xmax = tanh(upperCL)),  height = 0, show.legend = F, size = 1.2) +
  geom_vline(xintercept = 0, linetype = 2, colour = "black", alpha = 0.3) +
  # creating dots and different size (bee-swarm and bubbles)
  geom_point(aes(fill = name), size = 3, shape = 21) + #
  # setting colours
  scale_color_manual(values =   c("MutualistBoth"= colour_ls[1], "MutualistHorizontal"= colour_ls[2], "MutualistVertical" = colour_ls[3],"ParasiteBoth"= colour_ls[4], "ParasiteHorizontal" = colour_ls[5])) +
  scale_fill_manual(values = c("MutualistBoth"= colour_ls[1], "MutualistHorizontal"= colour_ls[2], "MutualistVertical" = colour_ls[3],"ParasiteBoth"= colour_ls[4], "ParasiteHorizontal" = colour_ls[5])) +
  scale_y_discrete(labels = c("MutualistBoth" = "Mutualist-\nBoth", "MutualistHorizontal" = "Mutualist-\nHorizontal","MutualistVertical" = "Mutualist-\nVertical", "ParasiteBoth" = "Parasite-\nBoth", "ParasiteHorizontal" = "Parasite-\nHorizontal")) +
  annotate('text', x = 0.93, y = 1:5 + 0.15, label= paste("italic(k)==", res_symbiosis_transmission1$n), parse=TRUE, hjust = "left", size=3.5) +
  labs(x = expression(paste(italic(r), " (correlation)")), y = "", size = expression(paste(italic(n), " (# of species)")) ,tag = "g" ) +
  guides(fill = "none", colour = "none") +
  theme_bw() +
  theme(legend.position="none") +
  theme(axis.text.y = element_text(size = 10, colour ="black",hjust = 0.5, angle = 90)) +
  # putting pictures in
  annotation_custom(rasterGrob(image_mutualism), xmin = -1.1, xmax = -0.9, ymin = 0.6, ymax = 1.2) + 
  annotation_custom(rasterGrob(image_both), xmin = -0.9, xmax = -0.6, ymin = 0.4, ymax = 1.2) + 
  annotation_custom(rasterGrob(image_mutualism), xmin = -1.1, xmax = -0.9, ymin = 1.6, ymax = 2.2) +
  annotation_custom(rasterGrob(image_horizontal),xmin = -0.9, xmax = -0.6, ymin = 1.4, ymax = 2.2) +
  annotation_custom(rasterGrob(image_mutualism), xmin = -1.1, xmax = -0.9, ymin = 2.6, ymax = 3.2) + 
  annotation_custom(rasterGrob(image_vertical), xmin = -0.9, xmax = -0.6, ymin = 2.4, ymax = 3.2) + 
  annotation_custom(rasterGrob(image_parasitism), xmin = -1.1, xmax = -0.9, ymin = 3.6, ymax = 4.2) +
  annotation_custom(rasterGrob(image_both), xmin = -0.9, xmax = -0.6, ymin = 3.4, ymax = 4.2) +
  annotation_custom(rasterGrob(image_parasitism), xmin = -1.1, xmax = -0.9, ymin = 4.6, ymax = 5.2) +
  annotation_custom(rasterGrob(image_horizontal), xmin = -0.9, xmax = -0.6, ymin = 4.4, ymax = 5.2)

Figure 3g: A forest plot showing the group-wise means (the categorical variable symbiosis_transmission) with their 95% confidences interval (thick lines) and 95% prediction intervals (thin lines) with observed effect sizes based on various sample sizes.

Putting together Figure 3

# building fig 3 using patchwork
fig3 <- (a/b/c/d + plot_layout(heights = c(1.6, 2, 3.7, 3.7))) | (e/f/g + plot_layout(heights = c(2.8, 
    2.8, 4.4)))  #+ plot_annotation(tag_levels = 'a', tag_suffix = ')')

fig3

# ggsave('../figs/fig3.png', width = 14, height = 14)
# ggsave('../figs/fig3.pdf', width = 14, height = 14)

Figure 3: putting all 7 panels together: Figure 3a - Figure 3g (see the main text)

Additional analyses (uni-predictor)

These are extra analyses not discussed in the main text.

The combined effect of host taxa and symbiosis (parasitism vs. mutualism)

# reordering
dat$host_tax_symbiosis <- factor(dat$host_tax_symbiosis, levels = c("MicrobeMutualist", 
    "MicrobeParasite", "PlantMutualist", "PlantParasite", "InvertMutualist", 
    "InvertParasite", "VertMutualist", "VertParasite"))

# meta-regression: mutiple intercepts
mr_host_tax_symbiosis1 <- rma.mv(yi = Zr, V = VZr, mods = ~host_tax_symbiosis - 
    1, test = "t", random = ~1 | authors, data = dat)

# # meta-regression: contrasts x 10 getting the level names out
level_names <- levels(dat$host_tax_symbiosis)

# helper function to run metafor meta-regression
run_rma <- function(name) {
    rma.mv(yi = Zr, V = VZr, mods = ~relevel(host_tax_symbiosis, ref = name), 
        test = "t", random = ~1 | authors, data = dat)
}

# results of meta-regression including all contrast results; taking the last
# level out ([-length(level_names)])
mr_host_tax_symbiosis <- map(level_names[-length(level_names)], run_rma)

Supplementary Table 10: Regression coefficients (Estimate), 95% confidence intervals (CIs), variance components (V) and variance explained, R²_[marginal] (R2) from the meta-regression with host_tax_symbiosis.

# getting marginal R2
r2_host_tax_symbiosis1 <- R2(mr_host_tax_symbiosis1)

# getting estimates
res_host_tax_symbiosis1 <- get_est(mr_host_tax_symbiosis1, mod = "host_tax_symbiosis")
res_host_tax_symbiosis <- map(mr_host_tax_symbiosis, ~get_est(.x, mod = "host_tax_symbiosis"))

# a list of the numbers to take out unnecessary contrasts
contra_list <- Map(seq, from = 1, to = 1:7)

# you need to flatten twice: first to make it a list and make it a vector
estimates <- map2(res_host_tax_symbiosis, contra_list, ~.x[-(.y), "estimate"]) %>% 
    flatten() %>% flatten_dbl()
lowerCLs <- map2(res_host_tax_symbiosis, contra_list, ~.x[-(.y), "lowerCL"]) %>% 
    flatten() %>% flatten_dbl()
upperCLs <- map2(res_host_tax_symbiosis, contra_list, ~.x[-(.y), "upperCL"]) %>% 
    flatten() %>% flatten_dbl()

# creating a table
tibble(`Fixed effect` = c(as.character(res_host_tax_symbiosis1$name), cont_gen(res_host_tax_symbiosis1$name)), 
    Estimate = c(res_host_tax_symbiosis1$estimate, estimates), `Lower CI [0.025]` = c(res_host_tax_symbiosis1$lowerCL, 
        lowerCLs), `Upper CI  [0.975]` = c(res_host_tax_symbiosis1$upperCL, 
        upperCLs), `V[authors]` = c(mr_host_tax_symbiosis1$sigma2, rep(NA, (8 + 
        choose(8, 2)) - 1)), R2 = c(r2_host_tax_symbiosis1[1], rep(NA, (8 + 
        choose(8, 2)) - 1))) %>% kable("html", digits = 3) %>% kable_styling("striped", 
    position = "left") %>% scroll_box(width = "100%", height = "300px")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	V[authors]	R2
MicrobeMutualist	0.806	0.596	1.015	0.027	0.296
MicrobeParasite	0.286	0.072	0.500	NA	NA
PlantMutualist	0.241	0.129	0.354	NA	NA
PlantParasite	0.246	0.155	0.337	NA	NA
InvertMutualist	0.424	0.348	0.500	NA	NA
InvertParasite	0.348	0.210	0.486	NA	NA
VertMutualist	0.674	0.143	1.206	NA	NA
VertParasite	0.321	0.272	0.371	NA	NA
MicrobeMutualist-MicrobeParasite	-0.520	-0.819	-0.220	NA	NA
MicrobeMutualist-PlantMutualist	-0.565	-0.802	-0.327	NA	NA
MicrobeMutualist-PlantParasite	-0.560	-0.788	-0.332	NA	NA
MicrobeMutualist-InvertMutualist	-0.382	-0.604	-0.159	NA	NA
MicrobeMutualist-InvertParasite	-0.458	-0.708	-0.207	NA	NA
MicrobeMutualist-VertMutualist	-0.131	-0.703	0.440	NA	NA
MicrobeMutualist-VertParasite	-0.484	-0.699	-0.269	NA	NA
MicrobeParasite-PlantMutualist	-0.045	-0.287	0.197	NA	NA
MicrobeParasite-PlantParasite	-0.040	-0.272	0.193	NA	NA
MicrobeParasite-InvertMutualist	0.138	-0.089	0.366	NA	NA
MicrobeParasite-InvertParasite	0.062	-0.193	0.317	NA	NA
MicrobeParasite-VertMutualist	0.389	-0.185	0.962	NA	NA
MicrobeParasite-VertParasite	0.036	-0.184	0.255	NA	NA
PlantMutualist-PlantParasite	0.005	-0.131	0.141	NA	NA
PlantMutualist-InvertMutualist	0.183	0.047	0.319	NA	NA
PlantMutualist-InvertParasite	0.107	-0.071	0.285	NA	NA
PlantMutualist-VertMutualist	0.433	-0.110	0.977	NA	NA
PlantMutualist-VertParasite	0.080	-0.042	0.203	NA	NA
PlantParasite-InvertMutualist	0.178	0.060	0.296	NA	NA
PlantParasite-InvertParasite	0.102	-0.063	0.267	NA	NA
PlantParasite-VertMutualist	0.428	-0.111	0.968	NA	NA
PlantParasite-VertParasite	0.075	-0.028	0.179	NA	NA
InvertMutualist-InvertParasite	-0.076	-0.234	0.082	NA	NA
InvertMutualist-VertMutualist	0.250	-0.287	0.788	NA	NA
InvertMutualist-VertParasite	-0.103	-0.193	-0.012	NA	NA
InvertParasite-VertMutualist	0.326	-0.223	0.876	NA	NA
InvertParasite-VertParasite	-0.027	-0.170	0.116	NA	NA
VertMutualist-VertParasite	-0.353	-0.887	0.181	NA	NA

# adding sample size (k) for each category
k_host_tax_symbiosis <- dat %>% group_by(host_tax_symbiosis) %>% count()
# getting estimates and predicitons
pred_host_tax_symbiosis <- get_pred(mr_host_tax_symbiosis1, mod = "host_tax_symbiosis") 
res_host_tax_symbiosis1 <- left_join(res_host_tax_symbiosis1, k_host_tax_symbiosis, by =  c("name" = "host_tax_symbiosis"))  %>% left_join(pred_host_tax_symbiosis)
#res_symbiosis1 
# drawing a funnel plot - fig 2b
fig_host_tax_symbiosis <- ggplot(data = res_host_tax_symbiosis1, aes(x = tanh(estimate), y = name)) +
  scale_x_continuous(limits=c(-1, 1), breaks = seq(-1, 1, by = 0.2) ) +
  geom_quasirandom(data = dat %>% filter(!is.na(host_tax_symbiosis)), 
                   aes(x= tanh(Zr), y = host_tax_symbiosis, size = ((1/VZr) + 3), colour = host_tax_symbiosis), groupOnX = FALSE, alpha=0.4) + 
  # 95 %precition interval (PI)
  geom_errorbarh(aes(xmin = tanh(lowerPR), xmax = tanh(upperPR)),  height = 0, show.legend = F, size = 0.5, alpha = 0.6) +
  # 95 %CI
  geom_errorbarh(aes(xmin = tanh(lowerCL), xmax = tanh(upperCL)),  height = 0, show.legend = F, size = 1.2) +
  geom_vline(xintercept = 0, linetype = 2, colour = "black", alpha = 0.3) +
  # creating dots and different size (bee-swarm and bubbles)
  geom_point(aes(fill = name), size = 3, shape = 21) + #
  # setting colours
  scale_color_manual(values =   c("MicrobeMutualist"= colour_ls[1], "MicrobeParasite"= colour_ls[2],  "PlantMutualist"= colour_ls[3], "PlantParasite"= colour_ls[4], "InvertMutualist" = colour_ls[5],  "InvertParasite"= colour_ls[6], "VertMutualist"= colour_ls[7],     "VertParasite"= colour_ls[8] )) +
  scale_fill_manual(values = c("MicrobeMutualist"= colour_ls[1], "MicrobeParasite"= colour_ls[2],  "PlantMutualist"= colour_ls[3], "PlantParasite"= colour_ls[4], "InvertMutualist" = colour_ls[5],  "InvertParasite"= colour_ls[6], "VertMutualist"= colour_ls[7],     "VertParasite"= colour_ls[8] )) +
  scale_y_discrete(labels = c("MicrobeMutualist"= "Microbe-\nMutualist", "MicrobeParasite"= "Microbe-\nParasite",  "PlantMutualist" = "Plant-\nMutualist", "PlantParasite"="Plant-\nParasite", "InvertMutualist" = "Invertebrate-\nMutualist",   "InvertParasite"= "Invertebrate-\nParasite", "VertMutualist"= "Vertebrate-\nMutualist",  "VertParasite"= "Vertebrate-\nParasite" )) +
  annotate('text', x = 0.93, y = 1:8 + 0.15, label= paste("italic(k)==", res_host_tax_symbiosis1$n), parse=TRUE, hjust = "left", size=3.5) +
  labs(x = expression(paste(italic(r), " (correlation)")), y = "", size = expression(paste(italic(n), " (# of species)")) ) +
  guides(fill = "none", colour = "none") +
  theme_bw() +
  theme(legend.position= c(0, 1), legend.justification = c(0,1)) +
  theme(legend.direction="horizontal") +
  #theme(legend.background = element_rect(fill = "white", colour = "black")) +
  theme(legend.background = element_blank()) +
  theme(axis.text.y = element_text(size = 10, colour ="black", hjust = 0.5, angle = 90)) +
  # putting pictures in
  annotation_custom(rasterGrob(image_microbe_host), xmin = -1.1, xmax = -0.9, ymin = 0.6, ymax = 1.2) + 
  annotation_custom(rasterGrob(image_mutualism), xmin = -0.9, xmax = -0.7, ymin = 0.6, ymax = 1.2) + 
  annotation_custom(rasterGrob(image_microbe_host), xmin = -1.1, xmax = -0.9, ymin = 1.6, ymax = 2.2) +
  annotation_custom(rasterGrob(image_parasitism),xmin = -0.9, xmax = -0.7, ymin = 1.6, ymax = 2.2) +
  annotation_custom(rasterGrob(image_plant_host), xmin = -1.1, xmax = -0.9, ymin = 2.6, ymax = 3.2) + 
  annotation_custom(rasterGrob(image_mutualism), xmin = -0.9, xmax = -0.7, ymin = 2.6, ymax = 3.2) + 
  annotation_custom(rasterGrob(image_plant_host), xmin = -1.1, xmax = -0.9, ymin = 3.6, ymax = 4.2) +
  annotation_custom(rasterGrob(image_parasitism), xmin = -0.9, xmax = -0.7, ymin = 3.6, ymax = 4.2) +
  annotation_custom(rasterGrob(image_invertebrate_host), xmin = -1.1, xmax = -0.9, ymin = 4.6, ymax = 5.2) +
  annotation_custom(rasterGrob(image_mutualism), xmin = -0.9, xmax = -0.7, ymin = 4.6, ymax = 5.2) +
  annotation_custom(rasterGrob(image_invertebrate_host), xmin = -1.1, xmax = -0.9, ymin = 5.6, ymax = 6.2) + 
  annotation_custom(rasterGrob(image_parasitism), xmin = -0.9, xmax = -0.7, ymin = 5.6, ymax = 6.2) + 
  annotation_custom(rasterGrob(image_vertebrate_host), xmin = -1.1, xmax = -0.9, ymin = 6.6, ymax = 7.2) +
  annotation_custom(rasterGrob(image_mutualism), xmin = -0.9, xmax = -0.7, ymin = 6.6, ymax = 7.2) +
  annotation_custom(rasterGrob(image_vertebrate_host), xmin = -1.1, xmax = -0.9, ymin = 7.6, ymax = 8.2) +
  annotation_custom(rasterGrob(image_parasitism), xmin = -0.9, xmax = -0.7, ymin = 7.6, ymax = 8.2)

fig_host_tax_symbiosis

Supplementary Figure 3: A forest plot showing the group-wise means (the categorical variable host_tax_symbiosis) with their 95% confidences interval (thick lines) and 95% prediction intervals (thin lines) with observed effect sizes based on various sample sizes.

Splitting host taxonomy by mode of symbiosis revealed that the observed higher phylogenetic congruence of host-symbiont cophylogenies involving a microbial host is driven primarily by greater congruence between microbial hosts and mutualist symbionts. Congruence is also relatively high for invertebrate hosts that harbour a mutualistic symbiont, while congruence appears to be lowest for plant hosts that harbour a parasitic symbiont.

The combined effect of symbiont taxa and symbiosis (parasitism vs. mutualism)

# reordering
dat$symbiont_tax_symbiosis <- factor(dat$symbiont_tax_symbiosis, levels = c("MicrobeMutualist", 
    "MicrobeParasite", "PlantMutualist", "PlantParasite", "InvertMutualist", 
    "InvertParasite", "VertParasite"))

# meta-regression: multiple intercepts
mr_symbiont_tax_symbiosis1 <- rma.mv(yi = Zr, V = VZr, mods = ~symbiont_tax_symbiosis - 
    1, test = "t", random = ~1 | authors, data = dat)

# # meta-regression: contrasts x 10 getting the level names out
level_names <- levels(dat$symbiont_tax_symbiosis)

# helper function to run metafor meta-regression
run_rma <- function(name) {
    rma.mv(yi = Zr, V = VZr, mods = ~relevel(symbiont_tax_symbiosis, ref = name), 
        test = "t", random = ~1 | authors, data = dat)
}

# results of meta-regression including all contrast results; taking the last
# level out ([-length(level_names)])
mr_symbiont_tax_symbiosis <- map(level_names[-length(level_names)], run_rma)

Supplementary Table 11: Regression coefficients (Estimate), 95% confidence intervals (CIs), variance components (V) and variance explained, R²_[marginal] (R2) from the meta-regression with symbiont_tax_symbiosis.

# getting marginal R2
r2_symbiont_tax_symbiosis1 <- R2(mr_symbiont_tax_symbiosis1)

# getting estimates
res_symbiont_tax_symbiosis1 <- get_est(mr_symbiont_tax_symbiosis1, mod = "symbiont_tax_symbiosis")
res_symbiont_tax_symbiosis <- map(mr_symbiont_tax_symbiosis, ~get_est(.x, mod = "symbiont_tax_symbiosis"))

# a list of the numbers to take out unnecessary contrasts
contra_list <- Map(seq, from = 1, to = 1:6)

# you need to flatten twice: first to make it a list and make it a vector
estimates <- map2(res_symbiont_tax_symbiosis, contra_list, ~.x[-(.y), "estimate"]) %>% 
    flatten() %>% flatten_dbl()
lowerCLs <- map2(res_symbiont_tax_symbiosis, contra_list, ~.x[-(.y), "lowerCL"]) %>% 
    flatten() %>% flatten_dbl()
upperCLs <- map2(res_symbiont_tax_symbiosis, contra_list, ~.x[-(.y), "upperCL"]) %>% 
    flatten() %>% flatten_dbl()

# creating a table
tibble(`Fixed effect` = c(as.character(res_symbiont_tax_symbiosis1$name), cont_gen(res_symbiont_tax_symbiosis1$name)), 
    Estimate = c(res_symbiont_tax_symbiosis1$estimate, estimates), `Lower CI [0.025]` = c(res_symbiont_tax_symbiosis1$lowerCL, 
        lowerCLs), `Upper CI  [0.975]` = c(res_symbiont_tax_symbiosis1$upperCL, 
        upperCLs), `V[authors]` = c(mr_symbiont_tax_symbiosis1$sigma2, rep(NA, 
        (7 + choose(7, 2)) - 1)), R2 = c(r2_symbiont_tax_symbiosis1[1], rep(NA, 
        (7 + choose(7, 2)) - 1))) %>% kable("html", digits = 3) %>% kable_styling("striped", 
    position = "left") %>% scroll_box(width = "100%", height = "300px")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	V[authors]	R2
MicrobeMutualist	0.437	0.365	0.509	0.029	0.263
MicrobeParasite	0.279	0.213	0.345	NA	NA
PlantMutualist	0.651	0.348	0.953	NA	NA
PlantParasite	1.499	0.555	2.443	NA	NA
InvertMutualist	0.290	0.165	0.415	NA	NA
InvertParasite	0.325	0.269	0.380	NA	NA
VertParasite	0.291	-0.100	0.681	NA	NA
MicrobeMutualist-MicrobeParasite	-0.158	-0.256	-0.061	NA	NA
MicrobeMutualist-PlantMutualist	0.214	-0.097	0.525	NA	NA
MicrobeMutualist-PlantParasite	1.062	0.115	2.008	NA	NA
MicrobeMutualist-InvertMutualist	-0.147	-0.291	-0.002	NA	NA
MicrobeMutualist-InvertParasite	-0.112	-0.203	-0.022	NA	NA
MicrobeMutualist-VertParasite	-0.146	-0.543	0.251	NA	NA
MicrobeParasite-PlantMutualist	0.372	0.062	0.682	NA	NA
MicrobeParasite-PlantParasite	1.220	0.274	2.166	NA	NA
MicrobeParasite-InvertMutualist	0.011	-0.130	0.153	NA	NA
MicrobeParasite-InvertParasite	0.046	-0.040	0.132	NA	NA
MicrobeParasite-VertParasite	0.012	-0.384	0.408	NA	NA
PlantMutualist-PlantParasite	0.848	-0.143	1.839	NA	NA
PlantMutualist-InvertMutualist	-0.361	-0.688	-0.033	NA	NA
PlantMutualist-InvertParasite	-0.326	-0.634	-0.019	NA	NA
PlantMutualist-VertParasite	-0.360	-0.854	0.134	NA	NA
PlantParasite-InvertMutualist	-1.208	-2.161	-0.256	NA	NA
PlantParasite-InvertParasite	-1.174	-2.120	-0.228	NA	NA
PlantParasite-VertParasite	-1.208	-2.229	-0.186	NA	NA
InvertMutualist-InvertParasite	0.034	-0.099	0.167	NA	NA
InvertMutualist-VertParasite	0.001	-0.409	0.411	NA	NA
InvertParasite-VertParasite	-0.034	-0.428	0.361	NA	NA

# adding sample size (k) for each category
k_symbiont_tax_symbiosis <- dat %>% group_by(symbiont_tax_symbiosis) %>% count()
# getting estimates and predicitons
pred_symbiont_tax_symbiosis <- get_pred(mr_symbiont_tax_symbiosis1, mod = "symbiont_tax_symbiosis") 
res_symbiont_tax_symbiosis1 <- left_join(res_symbiont_tax_symbiosis1, k_symbiont_tax_symbiosis, by =  c("name" = "symbiont_tax_symbiosis"))  %>% left_join(pred_symbiont_tax_symbiosis)
#res_symbiosis1 
# drawing a funnel plot - fig 2b
fig_symbiont_tax_symbiosis <- ggplot(data = res_symbiont_tax_symbiosis1, aes(x = tanh(estimate), y = name)) +
  scale_x_continuous(limits=c(-1, 1), breaks = seq(-1, 1, by = 0.2) ) +
  geom_quasirandom(data = dat %>% filter(!is.na(symbiont_tax_symbiosis)), 
                   aes(x= tanh(Zr), y = symbiont_tax_symbiosis, size = ((1/VZr) + 3), colour = symbiont_tax_symbiosis), groupOnX = FALSE, alpha=0.4) + 
  # 95 %precition interval (PI)
  geom_errorbarh(aes(xmin = tanh(lowerPR), xmax = tanh(upperPR)),  height = 0, show.legend = F, size = 0.5, alpha = 0.6) +
  # 95 %CI
  geom_errorbarh(aes(xmin = tanh(lowerCL), xmax = tanh(upperCL)),  height = 0, show.legend = F, size = 1.2) +
  geom_vline(xintercept = 0, linetype = 2, colour = "black", alpha = 0.3) +
  # creating dots and different size (bee-swarm and bubbles)
  geom_point(aes(fill = name), size = 3, shape = 21) + #
  # setting colours
  # setting colours
  scale_color_manual(values =   c("MicrobeMutualist"= colour_ls[1], "MicrobeParasite"= colour_ls[2],  "PlantMutualist"= colour_ls[3], "PlantParasite"= colour_ls[4], "InvertMutualist" = colour_ls[5],  "InvertParasite"= colour_ls[6], "VertParasite"= colour_ls[8] )) +
  scale_fill_manual(values = c("MicrobeMutualist"= colour_ls[1], "MicrobeParasite"= colour_ls[2],  "PlantMutualist"= colour_ls[3], "PlantParasite"= colour_ls[4], "InvertMutualist" = colour_ls[5],  "InvertParasite"= colour_ls[6], "VertParasite"= colour_ls[8] )) +
  scale_y_discrete(labels = c("MicrobeMutualist"= "Microbe-\nMutualist", "MicrobeParasite"= "Microbe-\nParasite",  "PlantMutualist" = "Plant-\nMutualist", "PlantParasite"="Plant-\nParasite", "InvertMutualist" = "Invertebrate-\nMutualist",   "InvertParasite"= "Invertebrate-\nParasite",  "VertParasite"= "Vertebrate-\nParasite" )) +
  annotate('text', x = 0.93, y = 1:7 + 0.15, label= paste("italic(k)==", res_symbiont_tax_symbiosis1$n), parse=TRUE, hjust = "left", size=3.5) +
  labs(x = expression(paste(italic(r), " (correlation)")), y = "", size = expression(paste(italic(n), " (# of species)")) ) +
  guides(fill = "none", colour = "none") +
  theme_bw() +
  theme(legend.position= c(0, 1), legend.justification = c(0, 1)) +
  theme(legend.direction="horizontal") +
  #theme(legend.background = element_rect(fill = "white", colour = "black")) +
  theme(legend.background = element_blank()) +
  theme(axis.text.y = element_text(size = 10, colour ="black", hjust = 0.5, angle = 90)) +
    # putting pictures in
  annotation_custom(rasterGrob(image_microbe_parasite), xmin = -1.1, xmax = -0.9, ymin = 0.6, ymax = 1.2) + 
  annotation_custom(rasterGrob(image_mutualism), xmin = -0.9, xmax = -0.7, ymin = 0.6, ymax = 1.2) + 
  annotation_custom(rasterGrob(image_microbe_parasite), xmin = -1.1, xmax = -0.9, ymin = 1.6, ymax = 2.2) +
  annotation_custom(rasterGrob(image_parasitism),xmin = -0.9, xmax = -0.7, ymin = 1.6, ymax = 2.2) +
  annotation_custom(rasterGrob(image_plant_parasite), xmin = -1.1, xmax = -0.9, ymin = 2.6, ymax = 3.2) + 
  annotation_custom(rasterGrob(image_mutualism), xmin = -0.9, xmax = -0.7, ymin = 2.6, ymax = 3.2) + 
  annotation_custom(rasterGrob(image_plant_parasite), xmin = -1.1, xmax = -0.9, ymin = 3.6, ymax = 4.2) +
  annotation_custom(rasterGrob(image_parasitism), xmin = -0.9, xmax = -0.7, ymin = 3.6, ymax = 4.2) +
  annotation_custom(rasterGrob(image_invertebrate_parasite), xmin = -1.1, xmax = -0.9, ymin = 4.6, ymax = 5.2) +
  annotation_custom(rasterGrob(image_mutualism), xmin = -0.9, xmax = -0.7, ymin = 4.6, ymax = 5.2) +
  annotation_custom(rasterGrob(image_invertebrate_parasite), xmin = -1.1, xmax = -0.9, ymin = 5.6, ymax = 6.2) + 
  annotation_custom(rasterGrob(image_parasitism), xmin = -0.9, xmax = -0.7, ymin = 5.6, ymax = 6.2) + 
  annotation_custom(rasterGrob(image_vertebrate_parasite), xmin = -1.1, xmax = -0.9, ymin = 6.6, ymax = 7.2) +
  annotation_custom(rasterGrob(image_parasitism), xmin = -0.9, xmax = -0.7, ymin = 6.6, ymax = 7.2) 

fig_symbiont_tax_symbiosis

Supplementary Figure 4: A forest plot showing the group-wise means (the categorical variable symbiont_tax_symbiosis) with their 95% confidences interval (thick lines) and 95% prediction intervals (thin lines) with observed effect sizes based on various sample sizes.

Splitting symbiont taxonomy by mode of symbiosis revealed much less variation, except for higher congruence exhibited by cophylogenies involving a plant symbiont (which are relatively rare), and the finding that cophylogenies involving a microbial mutualist symbiont are slightly more congruent than the remaining categories.

The combined effect of host and symbiont taxa

# reordering
dat$host_symbiont_tax <- factor(dat$host_symbiont_tax, levels = c("MicrobeInvert", 
    "MicrobeMicrobe", "MicrobePlant", "PlantInvert", "PlantMicrobe", "InvertInvert", 
    "InvertMicrobe", "InvertPlant", "VertInvert", "VertMicrobe", "VertVert"))

# meta-regression: multiple intercepts
mr_host_symbiont_tax1 <- rma.mv(yi = Zr, V = VZr, mods = ~host_symbiont_tax - 
    1, test = "t", random = ~1 | authors, data = dat)

# # meta-regression: contrasts x 10 getting the level names out
level_names <- levels(dat$host_symbiont_tax)

# helper function to run metafor meta-regression
run_rma <- function(name) {
    rma.mv(yi = Zr, V = VZr, mods = ~relevel(host_symbiont_tax, ref = name), 
        test = "t", random = ~1 | authors, data = dat)
}

# results of meta-regression including all contrast results; taking the last
# level out ([-length(level_names)])
mr_host_symbiont_tax <- map(level_names[-length(level_names)], run_rma)

Supplementary Table 12: Regression coefficients (Estimate), 95% confidence intervals (CIs), variance components (V) and variance explained, R²_[marginal] (R2) from the meta-regression with host_symbiont_tax.

# getting marginal R2
r2_host_symbiont_tax1 <- R2(mr_host_symbiont_tax1)

# getting estimates
res_host_symbiont_tax1 <- get_est(mr_host_symbiont_tax1, mod = "host_symbiont_tax")
res_host_symbiont_tax <- map(mr_host_symbiont_tax, ~get_est(.x, mod = "host_symbiont_tax"))

# a list of the numbers to take out unnecessary contrasts
contra_list <- Map(seq, from = 1, to = 1:10)

# you need to flatten twice: first to make it a list and make it a vector
estimates <- map2(res_host_symbiont_tax, contra_list, ~.x[-(.y), "estimate"]) %>% 
    flatten() %>% flatten_dbl()

lowerCLs <- map2(res_host_symbiont_tax, contra_list, ~.x[-(.y), "lowerCL"]) %>% 
    flatten() %>% flatten_dbl()

upperCLs <- map2(res_host_symbiont_tax, contra_list, ~.x[-(.y), "upperCL"]) %>% 
    flatten() %>% flatten_dbl()

# creating a table
tibble(`Fixed effect` = c(as.character(res_host_symbiont_tax1$name), cont_gen(res_host_symbiont_tax1$name)), 
    Estimate = c(res_host_symbiont_tax1$estimate, estimates), `Lower CI [0.025]` = c(res_host_symbiont_tax1$lowerCL, 
        lowerCLs), `Upper CI  [0.975]` = c(res_host_symbiont_tax1$upperCL, upperCLs), 
    `V[authors]` = c(mr_host_tax_symbiosis1$sigma2, rep(NA, (11 + choose(11, 
        2)) - 1)), R2 = c(r2_host_symbiont_tax1[1], rep(NA, (11 + choose(11, 
        2)) - 1))) %>% kable("html", digits = 3) %>% kable_styling("striped", 
    position = "left") %>% scroll_box(width = "100%", height = "300px")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	V[authors]	R2
MicrobeInvert	0.623	-0.119	1.366	0.027	0.22
MicrobeMicrobe	0.508	0.328	0.689	NA	NA
MicrobePlant	0.696	0.379	1.012	NA	NA
PlantInvert	0.216	0.119	0.312	NA	NA
PlantMicrobe	0.295	0.167	0.423	NA	NA
InvertInvert	0.421	0.286	0.556	NA	NA
InvertMicrobe	0.405	0.325	0.485	NA	NA
InvertPlant	0.929	0.187	1.672	NA	NA
VertInvert	0.355	0.287	0.422	NA	NA
VertMicrobe	0.293	0.215	0.371	NA	NA
VertVert	0.291	-0.107	0.688	NA	NA
MicrobeInvert-MicrobeMicrobe	-0.115	-0.879	0.649	NA	NA
MicrobeInvert-MicrobePlant	0.072	-0.735	0.880	NA	NA
MicrobeInvert-PlantInvert	-0.408	-1.157	0.341	NA	NA
MicrobeInvert-PlantMicrobe	-0.328	-1.082	0.425	NA	NA
MicrobeInvert-InvertInvert	-0.202	-0.957	0.553	NA	NA
MicrobeInvert-InvertMicrobe	-0.218	-0.965	0.529	NA	NA
MicrobeInvert-InvertPlant	0.306	-0.744	1.356	NA	NA
MicrobeInvert-VertInvert	-0.269	-1.014	0.477	NA	NA
MicrobeInvert-VertMicrobe	-0.330	-1.077	0.417	NA	NA
MicrobeInvert-VertVert	-0.332	-1.174	0.510	NA	NA
MicrobeMicrobe-MicrobePlant	0.187	-0.177	0.552	NA	NA
MicrobeMicrobe-PlantInvert	-0.293	-0.498	-0.088	NA	NA
MicrobeMicrobe-PlantMicrobe	-0.213	-0.435	0.008	NA	NA
MicrobeMicrobe-InvertInvert	-0.087	-0.313	0.139	NA	NA
MicrobeMicrobe-InvertMicrobe	-0.103	-0.301	0.095	NA	NA
MicrobeMicrobe-InvertPlant	0.421	-0.343	1.185	NA	NA
MicrobeMicrobe-VertInvert	-0.154	-0.347	0.039	NA	NA
MicrobeMicrobe-VertMicrobe	-0.215	-0.412	-0.018	NA	NA
MicrobeMicrobe-VertVert	-0.217	-0.654	0.219	NA	NA
MicrobePlant-PlantInvert	-0.480	-0.811	-0.149	NA	NA
MicrobePlant-PlantMicrobe	-0.401	-0.742	-0.059	NA	NA
MicrobePlant-InvertInvert	-0.274	-0.618	0.069	NA	NA
MicrobePlant-InvertMicrobe	-0.290	-0.617	0.036	NA	NA
MicrobePlant-InvertPlant	0.234	-0.574	1.041	NA	NA
MicrobePlant-VertInvert	-0.341	-0.665	-0.018	NA	NA
MicrobePlant-VertMicrobe	-0.403	-0.728	-0.077	NA	NA
MicrobePlant-VertVert	-0.405	-0.913	0.103	NA	NA
PlantInvert-PlantMicrobe	0.080	-0.081	0.240	NA	NA
PlantInvert-InvertInvert	0.206	0.040	0.372	NA	NA
PlantInvert-InvertMicrobe	0.190	0.064	0.315	NA	NA
PlantInvert-InvertPlant	0.714	-0.035	1.463	NA	NA
PlantInvert-VertInvert	0.139	0.021	0.257	NA	NA
PlantInvert-VertMicrobe	0.078	-0.047	0.202	NA	NA
PlantInvert-VertVert	0.075	-0.334	0.485	NA	NA
PlantMicrobe-InvertInvert	0.126	-0.060	0.312	NA	NA
PlantMicrobe-InvertMicrobe	0.110	-0.041	0.261	NA	NA
PlantMicrobe-InvertPlant	0.634	-0.119	1.388	NA	NA
PlantMicrobe-VertInvert	0.059	-0.085	0.204	NA	NA
PlantMicrobe-VertMicrobe	-0.002	-0.152	0.148	NA	NA
PlantMicrobe-VertVert	-0.004	-0.422	0.413	NA	NA
InvertInvert-InvertMicrobe	-0.016	-0.173	0.141	NA	NA
InvertInvert-InvertPlant	0.508	-0.247	1.263	NA	NA
InvertInvert-VertInvert	-0.067	-0.211	0.078	NA	NA
InvertInvert-VertMicrobe	-0.128	-0.284	0.028	NA	NA
InvertInvert-VertVert	-0.130	-0.550	0.290	NA	NA
InvertMicrobe-InvertPlant	0.524	-0.223	1.271	NA	NA
InvertMicrobe-VertInvert	-0.051	-0.156	0.054	NA	NA
InvertMicrobe-VertMicrobe	-0.112	-0.224	0.000	NA	NA
InvertMicrobe-VertVert	-0.114	-0.520	0.291	NA	NA
InvertPlant-VertInvert	-0.575	-1.320	0.171	NA	NA
InvertPlant-VertMicrobe	-0.636	-1.383	0.111	NA	NA
InvertPlant-VertVert	-0.638	-1.481	0.204	NA	NA
VertInvert-VertMicrobe	-0.061	-0.165	0.042	NA	NA
VertInvert-VertVert	-0.064	-0.467	0.340	NA	NA
VertMicrobe-VertVert	-0.002	-0.407	0.403	NA	NA

# colour list
#colour_ls <- c("#000000", "#E69F00", "#56B4E9", "#009E73",  "#F0E422",  "#0072B2",  "#D55E00", "#CC79A7", "#00008B", "#8B0A50", "#54FF9F", "#999999")

# adding sample size (k) for each category
k_host_symbiont_tax <- dat %>% group_by(host_symbiont_tax) %>% count()
# getting estimates and predicitons
pred_host_symbiont_tax <- get_pred(mr_host_symbiont_tax1, mod = "host_symbiont_tax") 
res_host_symbiont_tax1 <- left_join(res_host_symbiont_tax1, k_host_symbiont_tax, by =  c("name" = "host_symbiont_tax"))  %>% left_join(pred_host_symbiont_tax)
#res_symbiosis1 
# drawing a funnel plot - fig 2b
fig_host_symbiont_tax <- ggplot(data = res_host_symbiont_tax1, aes(x = tanh(estimate), y = name)) +
  scale_x_continuous(limits=c(-1, 1), breaks = seq(-1, 1, by = 0.2) ) +
  geom_quasirandom(data = dat %>% filter(!is.na(host_symbiont_tax)), 
                   aes(x= tanh(Zr), y = host_symbiont_tax, size = ((1/VZr) + 3), colour = host_symbiont_tax), groupOnX = FALSE, alpha=0.4) + 
  # 95 %precition interval (PI)
  geom_errorbarh(aes(xmin = tanh(lowerPR), xmax = tanh(upperPR)),  height = 0, show.legend = F, size = 0.5, alpha = 0.6) +
  # 95 %CI
  geom_errorbarh(aes(xmin = tanh(lowerCL), xmax = tanh(upperCL)),  height = 0, show.legend = F, size = 1.2) +
  geom_vline(xintercept = 0, linetype = 2, colour = "black", alpha = 0.3) +
  # creating dots and different size (bee-swarm and bubbles)
  geom_point(aes(fill = name), size = 3, shape = 21) + #
  # setting colours
  scale_color_manual(values =  c("MicrobeInvert" = colour_ls[1],  "MicrobeMicrobe"= colour_ls[2], "MicrobePlant" = colour_ls[3], "PlantInvert" = colour_ls[4],"PlantMicrobe" = colour_ls[5], "InvertInvert"  = colour_ls[6],  "InvertMicrobe" = colour_ls[7], "InvertPlant" = colour_ls[8],"VertInvert"  = colour_ls[9], "VertMicrobe"= colour_ls[10],"VertVert"  = colour_ls[11])) +
  scale_fill_manual(values = c("MicrobeInvert" = colour_ls[1],  "MicrobeMicrobe"= colour_ls[2], "MicrobePlant" = colour_ls[3], "PlantInvert" = colour_ls[4],"PlantMicrobe" = colour_ls[5], "InvertInvert"  = colour_ls[6],  "InvertMicrobe" = colour_ls[7], "InvertPlant" = colour_ls[8],"VertInvert"  = colour_ls[9], "VertMicrobe"= colour_ls[10],"VertVert"  = colour_ls[11])) +
  scale_y_discrete(labels = c("MicrobeInvert" = "Microbe-\nInvertebrate",  "MicrobeMicrobe"= "Microbe-\nMicrobe", "MicrobePlant" = "Microbe-\nPlant", "PlantInvert" = "Plant-\nInvertebrate","PlantMicrobe" = "Plant-\nMicrobe", "InvertInvert"  = "Invertebrate\nInvertebrate",  "InvertMicrobe" = "Invertebrate-\nMicrobe", "InvertPlant" = "Invertebrate-\nPlant","VertInvert"  = "Vertebrate-\nInvertebrate", "VertMicrobe"= "Vertebrate-\nMicrobe", "VertVert"  = "Vertebrate-\nVertebrate")) +
  annotate('text', x = 0.93, y = 1:11 + 0.15, label= paste("italic(k)==", res_host_symbiont_tax1$n), parse=TRUE, hjust = "left", size=3.5) +
  labs(x = expression(paste(italic(r), " (correlation)")), y = "", size = expression(paste(italic(n), " (# of species)")) ) +
  guides(fill = "none", colour = "none") +
  theme_bw() +
  theme(legend.position= c(0, 1), legend.justification = c(0,1)) +
  theme(legend.direction="horizontal") +
  #theme(legend.background = element_rect(fill = "white", colour = "black")) +
  theme(legend.background = element_blank()) +
  theme(axis.text.y = element_text(size = 10, colour ="black", hjust = 0.5, angle = 90)) +
    # putting pictures in
  annotation_custom(rasterGrob(image_microbe_host), xmin = -1.1, xmax = -0.9, ymin = 0.6, ymax = 1.2) + 
  annotation_custom(rasterGrob(image_invertebrate_parasite), xmin = -0.9, xmax = -0.7, ymin = 0.6, ymax = 1.2) + 
  annotation_custom(rasterGrob(image_microbe_host), xmin = -1.1, xmax = -0.9, ymin = 1.6, ymax = 2.2) +
  annotation_custom(rasterGrob(image_microbe_parasite),xmin = -0.9, xmax = -0.7, ymin = 1.6, ymax = 2.2) +
  annotation_custom(rasterGrob(image_microbe_host), xmin = -1.1, xmax = -0.9, ymin = 2.6, ymax = 3.2) + 
  annotation_custom(rasterGrob(image_plant_parasite), xmin = -0.9, xmax = -0.7, ymin = 2.6, ymax = 3.2) + 
  #
  annotation_custom(rasterGrob(image_plant_host), xmin = -1.1, xmax = -0.9, ymin = 3.6, ymax = 4.2) +
  annotation_custom(rasterGrob(image_invertebrate_parasite), xmin = -0.9, xmax = -0.7, ymin = 3.6, ymax = 4.2) +
  annotation_custom(rasterGrob(image_plant_host), xmin = -1.1, xmax = -0.9, ymin = 4.6, ymax = 5.2) +
  annotation_custom(rasterGrob(image_microbe_parasite), xmin = -0.9, xmax = -0.7, ymin = 4.6, ymax = 5.2) +
  #
  annotation_custom(rasterGrob(image_invertebrate_host), xmin = -1.1, xmax = -0.9, ymin = 5.6, ymax = 6.2) + 
  annotation_custom(rasterGrob(image_invertebrate_parasite), xmin = -0.9, xmax = -0.7, ymin = 5.6, ymax = 6.2) + 
  annotation_custom(rasterGrob(image_invertebrate_host), xmin = -1.1, xmax = -0.9, ymin = 6.6, ymax = 7.2) +
  annotation_custom(rasterGrob(image_microbe_parasite), xmin = -0.9, xmax = -0.7, ymin = 6.6, ymax = 7.2) +
  annotation_custom(rasterGrob(image_invertebrate_host), xmin = -1.1, xmax = -0.9, ymin = 7.6, ymax = 8.2) +
  annotation_custom(rasterGrob(image_plant_parasite), xmin = -0.9, xmax = -0.7, ymin = 7.6, ymax = 8.2) +
  #
  annotation_custom(rasterGrob(image_vertebrate_host), xmin = -1.1, xmax = -0.9, ymin = 8.6, ymax = 9.2) + 
  annotation_custom(rasterGrob(image_invertebrate_parasite), xmin = -0.9, xmax = -0.7, ymin = 8.6, ymax = 9.2) + 
  annotation_custom(rasterGrob(image_vertebrate_host), xmin = -1.1, xmax = -0.9, ymin = 9.6, ymax = 10.2) +
  annotation_custom(rasterGrob(image_microbe_parasite), xmin = -0.9, xmax = -0.7, ymin = 9.6, ymax = 10.2) +
  annotation_custom(rasterGrob(image_vertebrate_host), xmin = -1.1, xmax = -0.9, ymin = 10.6, ymax = 11.2) +
  annotation_custom(rasterGrob(image_vertebrate_parasite), xmin = -0.9, xmax = -0.7, ymin = 10.6, ymax = 11.2)


fig_host_symbiont_tax

Supplementary Figure 4: A forest plot showing the group-wise means (the categorical variable host_symbiont_tax) with their 95% confidences interval (thick lines) and 95% prediction intervals (thin lines) with observed effect sizes based on various sample sizes.

Model selection (multi-predictor model)

Here we build the best model via an AICc based model selection method implemented in the R package MuMin¹⁰. For the full model, we had 6 variables: symbiosis, host_tax_broad, symbiont_tax_broad, mode_of_transmission_broad, endo_or_ecto, & log(host_range_link_ratio). We did not use log(host_range_taxonomic_breadth) as it is co-linear with log(host_range_link_ratio) and also any of interaction terms.

# creates a new function to run in MuMIn
updated.rma.mv <- updateable(rma.mv)
# updated.rma.mv

# testing the new function use method = 'ML' so that we can compare AIC
mr_full <- updated.rma.mv(yi = Zr, V = VZr, mods = ~symbiosis + host_tax_broad + 
    symbiont_tax_broad + mode_of_transmission_broad + endo_or_ecto + log(host_range_link_ratio), 
    test = "t", random = ~1 | authors, method = "ML", data = dat)

# ============================= additional methods for 'rma.mv' class (made
# by Kamil Barton) we need this to run model selection with rma.mv in MuMIn
# =============================
formula.rma.mv <- function(x, ...) return(eval(getCall(x)$mods))

makeArgs.rma.mv <- function(obj, termNames, comb, opt, ...) {
    ret <- MuMIn:::makeArgs.default(obj, termNames, comb, opt)
    names(ret)[1L] <- "mods"
    ret
}

nobs.rma.mv <- function(object, ...) attr(logLik(object), "nall")

coefTable.rma.mv <- function(model, ...) MuMIn:::.makeCoefTable(model$b, model$se, 
    coefNames = rownames(model$b))
# =============================

# testing dredge dredge(full.model, evaluate=F) # show all candidate models
# n = 32 model exisit
candidates <- dredge(mr_full)

# displays delta AICc <2
candidates_aic2 <- subset(candidates, delta < 2)

# model averaging it seems like models are using z values rather than t
# values (which will be OK)
mr_averaged_aic2 <- summary(model.avg(candidates, delta < 2))

# relative importance of each predictor
importance <- importance(candidates)

# use REML if not for model comparision
model1 <- rma.mv(yi = Zr, V = VZr, mods = ~mode_of_transmission_broad + host_tax_broad, 
    test = "t", random = ~1 | authors, method = "REML", data = dat)
model2 <- rma.mv(yi = Zr, V = VZr, mods = ~mode_of_transmission_broad + host_tax_broad + 
    symbiosis, test = "t", random = ~1 | authors, method = "REML", data = dat)

Supplementary Table 13: The top 2 models (out of 32 possible models) within the \(\Delta\)AIC difference of 2 and which 6 variables: symbiosis, host_tax_broad, symbiont_tax_broad, mode_of_transmission_broad, endo_or_ecto, & log(host_range_link_ratio) were included (indicated by \(+\)); model weights (for the 2 models) and the sum of weights for each of the variable (from the 32 models) are included.

# creating a table
tibble(`Model (variable weight)` = c("Model1", "Model2", "(Sum of weights)"), 
    transmission = c(if_else(candidates_aic2$mode_of_transmission_broad == "+", 
        "$+$", "NA"), round(importance[1], 3)), host_tax = c(if_else(candidates_aic2$host_tax_broad == 
        "+", "$+$", "NA"), round(importance[2], 3)), symbiosis = c(if_else(candidates_aic2$symbiosis == 
        "+", "$+$", "NA"), round(importance[3], 3)), symbiont_tax = c(if_else(candidates_aic2$symbiont_tax_broad == 
        "+", "$+$", "NA"), round(importance[4], 3)), endo_or_ecto = c(if_else(candidates_aic2$endo_or_ecto == 
        "+", "$+$", "NA"), round(importance[5], 3)), host_range = c(if_else(candidates_aic2$`log(host_range_link_ratio)` == 
        "+", "$+$", "NA"), round(importance[6], 3)), delta_AICc = c(candidates_aic2$delta, 
        NA), Weight = c(candidates_aic2$weight, NA)) %>% kable("html", digits = 3) %>% 
    kable_styling("striped", position = "left")

Model (variable weight)	transmission	host_tax	symbiosis	symbiont_tax	endo_or_ecto	host_range	delta_AICc	Weight
Model1	\(+\)	\(+\)	NA	NA	NA	NA	0.000	0.633
Model2	\(+\)	\(+\)	\(+\)	NA	NA	NA	1.087	0.367
(Sum of weights)	0.998	0.682	0.342	0.217	0.129	0.087	NA	NA

Model averaging

Supplementary Table 14: The averaged estimates for regression coefficients (Estimate), 95% confidence intervals (CIs), variance components (V) and variance explained, R²_[marginal] (R2) from the 2 best meta-regression models.

# getting averaged R2 and variance components not provided by the MuMIn
# package
average_sigma2 <- weighted.mean(x = c(model1$sigma2, model2$sigma2), w = candidates_aic2$weight)
average_R2 <- weighted.mean(x = c(R2(model1)[1], R2(model2)[1]), w = candidates_aic2$weight)

# creating a table
tibble(`Fixed effect` = c("Intercept (both-Microbe-Mutulist)", "Microbe-Plant", 
    "Microbe-Invert", "Microbe-Vert", "both-horizontal", "both-vertical", "Mutulist-Parasite"), 
    Estimate = mr_averaged_aic2$coefmat.full[, 1], `Lower CI [0.025]` = mr_averaged_aic2$coefmat.full[, 
        1] - mr_averaged_aic2$coefmat.full[, 2] * qnorm(0.975), `Upper CI  [0.975]` = mr_averaged_aic2$coefmat.full[, 
        1] + mr_averaged_aic2$coefmat.full[, 2] * qnorm(0.975), `V[authors]` = c(average_sigma2, 
        rep(NA, 6)), R2 = c(average_R2, rep(NA, 6))) %>% kable("html", digits = 3) %>% 
    kable_styling("striped", position = "left")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	V[authors]	R2
Intercept (both-Microbe-Mutulist)	0.547	0.382	0.711	0.024	0.282
Microbe-Plant	-0.248	-0.417	-0.078	NA	NA
Microbe-Invert	-0.185	-0.355	-0.015	NA	NA
Microbe-Vert	-0.173	-0.332	-0.014	NA	NA
both-horizontal	-0.055	-0.141	0.032	NA	NA
both-vertical	0.120	-0.031	0.271	NA	NA
Mutulist-Parasite	-0.016	-0.095	0.062	NA	NA

Publication Bias Analysis

Here, we conducted 3 kinds of publication bias analyses: 1) contour-enhanced funnel plots¹¹ of residuals^7,12, 2) a type of Egger regression^12,13, and 3) a regression-based time-lag bias test⁷.

Funnel plot

A normal funnel plot assumes the homogeneity (i.e., I² = 0). therefore, we controlled for important moderators (i.e., mode_of_transmission_broad, host_tax_broad, & symbiosis).

Residual funnel plot 1

It is not normal skewness that we see in our enhanced-counter funnel plot (Supplementary Figure 5). This funnel asymmetry seems different from one cased by publication bias¹¹; we do not expect a “hollow” in the region with high precision or i.e. inverse standard error (4.5-6) and relative high effect sizes, (Zr = 0.2-0.7). The funnel asymmetry is mainly caused by the boundary created (see “Sensitivity Analysis” section where we deal with this skewness).

# 
res_funnel_plot <- rma.mv(yi = Zr, V = VZr, mods = ~mode_of_transmission_broad + 
    host_tax_broad + symbiosis, random = ~1 | authors, data = dat)

funnel(res_funnel_plot, yaxis = "seinv", level = c(90, 95, 99), shade = c("white", 
    "gray55", "gray75"), refline = 0, legend = TRUE)

Supplementary Figure 5: A residual funnel plot from the meta-regression model with mode_of_transmission_broad, host_tax_broad, & symbiosis; ‘residual value’ is on Zr and ‘inverse standard error’ is precision 1/sqrt(VZr).

Residual funnel plot 2

Further, Egger regression analyses (see below) showed sqrt(VZr) (sampling errors [SE] for effect sizes) accounts much heterogeneity so that we added to our model. The funnel asymmetry we see in Supplementary Figure 6 (if any) is much less sever than that in Supplementary Figure 5.

# 
res_funnel_plot2 <- rma.mv(yi = Zr, V = VZr, mods = ~sqrt(VZr) + mode_of_transmission_broad + 
    host_tax_broad + symbiosis, random = ~1 | authors, data = dat)

funnel(res_funnel_plot2, yaxis = "seinv", level = c(90, 95, 99), shade = c("white", 
    "gray55", "gray75"), refline = 0, legend = TRUE)

Supplementary Figure 6: A residual funnel plot from the meta-regression model with sqrt(VZr), mode_of_transmission_broad, host_tax_broad, & symbiosis; ‘residual value’ is on Zr and ‘inverse standard error’ is precision 1/sqrt(VZr).

Egger regression

Egger regression tests whether funnel asymmetries we see in funnel plots are statistical significant or not.

Univariate Egger regression

The test (or sqrt(VZr)) is significant. However, as mentioned above, this is due to the boundary created by the number of randamizations; this boundary can be seen in Supplementary Figure 7

# 
egger_regression_uni <- rma.mv(yi = Zr, V = VZr, mods = ~sqrt(VZr), random = ~1 | 
    authors, data = dat)

Supplementary Table 15: Regression coefficients (Estimate), 95% confidence intervals (CIs), variance components (V) and variance explained, R²_[marginal] (R2) from the meta-regression with symbiont_tax_broad.

# getting marginal R2
r2_egger_regression_uni <- R2(egger_regression_uni)

# getting estimates: name does not work for slopes
res_egger_regression_uni <- get_est(egger_regression_uni, mod = "sqrt(VZr)")

# creating a table
tibble(`Fixed effect` = c("Intercept", "sqrt(VZr)"), Estimate = c(res_egger_regression_uni$estimate), 
    `Lower CI [0.025]` = c(res_egger_regression_uni$lowerCL), `Upper CI  [0.975]` = c(res_egger_regression_uni$upperCL), 
    `V[authors]` = c(egger_regression_uni$sigma2, NA), R2 = c(r2_egger_regression_uni[1], 
        NA)) %>% kable("html", digits = 3) %>% kable_styling("striped", position = "left")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	V[authors]	R2
Intercept	0.124	0.035	0.212	0.026	0.255
sqrt(VZr)	1.314	0.817	1.810	NA	NA

pred_egger_regression_uni <- predict.rma(egger_regression_uni)

# plotting

fit_egger_regression_uni <- dat %>% mutate(ymin = pred_egger_regression_uni$ci.lb, 
    ymax = pred_egger_regression_uni$ci.ub, ymin2 = pred_egger_regression_uni$cr.lb, 
    ymax2 = pred_egger_regression_uni$cr.ub, pred = pred_egger_regression_uni$pred) %>% 
    ggplot(aes(x = sqrt(VZr), y = Zr, size = (1/VZr) + 3)) + geom_point(shape = 21, 
    fill = "grey90") + # geom_ribbon(aes(ymin = ymin, ymax = ymax), fill = '#0072B2') + # not quite
# sure why this does not work
geom_smooth(aes(y = ymin2), method = "loess", se = FALSE, lty = "dotted", lwd = 0.25, 
    colour = "#0072B2") + geom_smooth(aes(y = ymax2), method = "loess", se = FALSE, 
    lty = "dotted", lwd = 0.25, colour = "#0072B2") + geom_smooth(aes(y = ymin), 
    method = "loess", se = FALSE, lty = "dotted", lwd = 0.25, colour = "#D55E00") + 
    geom_smooth(aes(y = ymax), method = "loess", se = FALSE, lty = "dotted", 
        lwd = 0.25, colour = "#D55E00") + geom_smooth(aes(y = pred), method = "loess", 
    se = FALSE, lty = "dashed", lwd = 0.5, colour = "black") + ylim(-1, 2) + 
    xlim(0.05, 0.45) + # geom_abline(intercept = mr_host_range_link_ratio$beta[[1]], slope =
# mr_host_range_link_ratio$beta[[2]], alpha = 0.7, linetype = 'dashed', size
# = 0.5) +
labs(x = "sqrt(sampling variance)", y = expression(paste(italic(Zr), " (effect size)")), 
    size = expression(paste(italic(n), " (# of species)"))) + guides(fill = "none", 
    colour = "none") + # themses
theme_bw() + theme(legend.position = c(0, 1), legend.justification = c(0, 1)) + 
    theme(legend.direction = "horizontal") + # theme(legend.background = element_rect(fill = 'white', colour = 'black'))
# +
theme(legend.background = element_blank()) + theme(axis.text.y = element_text(size = 10, 
    colour = "black", hjust = 0.5, angle = 90))

fit_egger_regression_uni

Supplementary Figure 7: A bubble plot showing a predicted regression line for the contentious variable log(host_range_link_ratio) with their 95% confidences regions (orange dotted lines) and 95% prediction regions (blue dotted lines) with observed effect sizes based on various sample sizes.

Multivariate Egger regression

We also conducted a Egger regression controlling other important moderators (i.e., mode_of_transmission_broad, host_tax_broad, & symbiosis). After controlling for these variables, sqrt(VZr) stays significant.

# 
egger_regression_mul <- rma.mv(yi = Zr, V = VZr, mods = ~sqrt(VZr) + mode_of_transmission_broad + 
    host_tax_broad + symbiosis, random = ~1 | authors, data = dat)

Supplementary Table 16: Regression coefficients (Estimate), 95% confidence intervals (CIs), variance components (V) and variance explained, R²_[marginal] (R2) from the meta-regression with symbiont_tax_broad.

# getting marginal R2
r2_egger_regression_mul <- R2(egger_regression_mul)

# creating a table
tibble(`Fixed effect` = c("Intercept (both-Microbe-Mutulist)", "sqrt(VZr)", 
    "both-horizontal", "both-vertical", "Microbe-Plant", "Microbe-Invert", "Microbe-Vert", 
    "Mutulist-Parasite"), Estimate = c(egger_regression_mul$b), `Lower CI [0.025]` = c(egger_regression_mul$ci.lb), 
    `Upper CI  [0.975]` = c(egger_regression_mul$ci.ub), `V[authors]` = c(egger_regression_mul$sigma2, 
        rep(NA, 7)), R2 = c(r2_egger_regression_mul[1], rep(NA, 7))) %>% kable("html", 
    digits = 3) %>% kable_styling("striped", position = "left")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	V[authors]	R2
Intercept (both-Microbe-Mutulist)	0.290	0.094	0.486	0.018	0.478
sqrt(VZr)	1.202	0.726	1.678	NA	NA
both-horizontal	-0.043	-0.125	0.039	NA	NA
both-vertical	0.138	-0.016	0.291	NA	NA
Microbe-Plant	-0.167	-0.331	-0.002	NA	NA
Microbe-Invert	-0.145	-0.309	0.019	NA	NA
Microbe-Vert	-0.095	-0.251	0.062	NA	NA
Mutulist-Parasite	-0.040	-0.144	0.065	NA	NA

pred_egger_regression_mul <-predict.rma(egger_regression_mul) 

# plotting

fit_egger_regression_mul <-  dat %>% 
  filter(!is.na(mode_of_transmission_broad) & !is.na(host_tax_broad) & !is.na(symbiosis))  %>% # getting ride of NA values
  mutate(ymin = pred_egger_regression_mul$ci.lb, 
         ymax = pred_egger_regression_mul$ci.ub,
         ymin2 = pred_egger_regression_mul$cr.lb,
         ymax2 = pred_egger_regression_mul$cr.ub,
         pred = pred_egger_regression_mul$pred) %>% 
  ggplot(aes(x = sqrt(VZr), y = Zr, size = (1/VZr) + 3)) +
  geom_point(shape = 21, fill = "grey90") +
  #geom_ribbon(aes(ymin = ymin, ymax = ymax), fill = "#0072B2")  + # not quite sure why this does not work
  geom_smooth(aes(y = ymin2), method =  "loess", se = FALSE, lty =  "dotted", lwd = 0.25, colour = "#0072B2") +
  geom_smooth(aes(y = ymax2), method =  "loess", se = FALSE, lty = "dotted", lwd = 0.25, colour = "#0072B2") +
  geom_smooth(aes(y = ymin), method =  "loess", se = FALSE,lty = "dotted", lwd = 0.25, colour ="#D55E00") +
  geom_smooth(aes(y = ymax), method =  "loess", se = FALSE, lty ="dotted", lwd = 0.25, colour ="#D55E00") + 
  geom_smooth(aes(y = pred), method =  "loess", se = FALSE, lty ="dashed", lwd = 0.5, colour ="black") +  
  ylim(-1, 2) + xlim(0.05, 0.45) +
  #geom_abline(intercept = mr_host_range_link_ratio$beta[[1]], slope = mr_host_range_link_ratio$beta[[2]], alpha = 0.7, linetype = "dashed", size = 0.5) +
  labs(x = "sqrt(sampling variance)", y = expression(paste(italic(Zr), " (effect size)")), size = expression(paste(italic(n), " (# of species)"))) +
  guides(fill = "none", colour = "none") +
  # themses
  theme_bw() +
  theme(legend.position= c(0, 1), legend.justification = c(0, 1)) +
  theme(legend.direction="horizontal") +
  #theme(legend.background = element_rect(fill = "white", colour = "black")) +
  theme(legend.background = element_blank()) +
  theme(axis.text.y = element_text(size = 10, colour ="black", hjust = 0.5, angle = 90)) 

fit_egger_regression_mul

Supplementary Figure 8: A bubble plot showing a predicted loess line for the contentious variable log(host_range_link_ratio) (given the values of the other 3 variables in the model) with their 95% confidences regions (orange dotted lines) and 95% prediction regions (blue dotted lines) with observed effect sizes based on various sample sizes. Note that the lines are not linear as these are based on multivariate predictions of the data points.

Time-lag bias

We do not find any evidence of the time-lag effect (a decline in the magnitude of the effect over time) in both the univariate and multivariate models (Supplementary Figure 9 and Supplementary Figure 10).

Univariate time-lag bias

# 
time_lag_effect_uni <- rma.mv(yi = Zr, V = VZr, mods = ~year, random = ~1 | 
    authors, data = dat)

Supplementary Table 17: Regression coefficients (Estimate), 95% confidence intervals (CIs), variance components (V) and variance explained, R²_[marginal] (R2) from the meta-regression with year.

# getting marginal R2
r2_time_lag_effect_uni <- R2(time_lag_effect_uni)

# getting estimates: name does not work for slopes
res_time_lag_effect_uni <- get_est(time_lag_effect_uni, mod = "year")

# creating a table
tibble(`Fixed effect` = c("Intercept", "Year"), Estimate = c(res_time_lag_effect_uni$estimate), 
    `Lower CI [0.025]` = c(res_time_lag_effect_uni$lowerCL), `Upper CI  [0.975]` = c(res_time_lag_effect_uni$upperCL), 
    `V[authors]` = c(time_lag_effect_uni$sigma2, NA), R2 = c(r2_time_lag_effect_uni[1], 
        NA)) %>% kable("html", digits = 3) %>% kable_styling("striped", position = "left")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	V[authors]	R2
Intercept	-2.547	-14.985	9.892	0.033	0.002
Year	0.001	-0.005	0.008	NA	NA

pred_time_lag_effect_uni <- predict.rma(time_lag_effect_uni)

# plotting

fit_time_lag_effect <- dat %>% mutate(ymin = pred_time_lag_effect_uni$ci.lb, 
    ymax = pred_time_lag_effect_uni$ci.ub, ymin2 = pred_time_lag_effect_uni$cr.lb, 
    ymax2 = pred_time_lag_effect_uni$cr.ub, pred = pred_time_lag_effect_uni$pred) %>% 
    ggplot(aes(x = year, y = Zr, size = (1/VZr) + 3)) + geom_point(shape = 21, 
    fill = "grey90") + # geom_ribbon(aes(ymin = ymin, ymax = ymax), fill = '#0072B2') + # not quite
# sure why this does not work
geom_smooth(aes(y = ymin2), method = "loess", se = FALSE, lty = "dotted", lwd = 0.25, 
    colour = "#0072B2") + geom_smooth(aes(y = ymax2), method = "loess", se = FALSE, 
    lty = "dotted", lwd = 0.25, colour = "#0072B2") + geom_smooth(aes(y = ymin), 
    method = "loess", se = FALSE, lty = "dotted", lwd = 0.25, colour = "#D55E00") + 
    geom_smooth(aes(y = ymax), method = "loess", se = FALSE, lty = "dotted", 
        lwd = 0.25, colour = "#D55E00") + geom_smooth(aes(y = pred), method = "loess", 
    se = FALSE, lty = "dashed", lwd = 0.5, colour = "black") + ylim(-1, 2) + 
    xlim(1994, 2019) + scale_x_continuous(breaks = c(1995, 2000, 2005, 2010, 
    2015, 2020)) + # geom_abline(intercept = mr_host_range_link_ratio$beta[[1]], slope =
# mr_host_range_link_ratio$beta[[2]], alpha = 0.7, linetype = 'dashed', size
# = 0.5) +
labs(x = "Year", y = expression(paste(italic(Zr), " (effect size)")), size = expression(paste(italic(n), 
    " (# of species)"))) + guides(fill = "none", colour = "none") + # themses
theme_bw() + theme(legend.position = c(0, 1), legend.justification = c(0, 1)) + 
    theme(legend.direction = "horizontal") + # theme(legend.background = element_rect(fill = 'white', colour = 'black'))
# +
theme(legend.background = element_blank()) + theme(axis.text.y = element_text(size = 10, 
    colour = "black", hjust = 0.5, angle = 90))

fit_time_lag_effect

Supplementary Figure 9: A bubble plot showing a predicted regression line for the contentious variable year with their 95% confidences regions (orange dotted lines) and 95% prediction regions (blue dotted lines) with observed effect sizes based on various sample sizes.

Multivariate time-lag bias

# 
time_lag_effect_mul <- rma.mv(yi = Zr, V = VZr, mods = ~year + mode_of_transmission_broad + 
    host_tax_broad + symbiosis, random = ~1 | authors, data = dat)

Supplementary Table 18: Regression coefficients (Estimate), 95% confidence intervals (CIs), variance components (V) and variance explained, R²_[marginal] (R2) from the meta-regression with year.

# getting marginal R2
r2_time_lag_effect_mul <- R2(time_lag_effect_mul)

# creating a table
tibble(`Fixed effect` = c("Intercept (both-Microbe-Mutulist)", "Year", "both-horizontal", 
    "both-vertical", "Microbe-Plant", "Microbe-Invert", "Microbe-Vert", "Mutulist-Parasite"), 
    Estimate = c(time_lag_effect_mul$b), `Lower CI [0.025]` = c(time_lag_effect_mul$ci.lb), 
    `Upper CI  [0.975]` = c(time_lag_effect_mul$ci.ub), `V[authors]` = c(time_lag_effect_mul$sigma2, 
        rep(NA, 7)), R2 = c(r2_time_lag_effect_mul[1], rep(NA, 7))) %>% kable("html", 
    digits = 3) %>% kable_styling("striped", position = "left")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	V[authors]	R2
Intercept (both-Microbe-Mutulist)	1.221	-10.447	12.890	0.024	0.28
Year	0.000	-0.006	0.005	NA	NA
both-horizontal	-0.053	-0.141	0.036	NA	NA
both-vertical	0.101	-0.064	0.265	NA	NA
Microbe-Plant	-0.255	-0.428	-0.082	NA	NA
Microbe-Invert	-0.187	-0.361	-0.014	NA	NA
Microbe-Vert	-0.168	-0.333	-0.003	NA	NA
Mutulist-Parasite	-0.045	-0.156	0.067	NA	NA

pred_time_lag_effect_mul <- predict.rma(time_lag_effect_mul)

# plotting
fit_time_lag_effect_mul <- dat %>% filter(!is.na(mode_of_transmission_broad) & 
    !is.na(host_tax_broad) & !is.na(symbiosis)) %>% mutate(ymin = pred_time_lag_effect_mul$ci.lb, 
    ymax = pred_time_lag_effect_mul$ci.ub, ymin2 = pred_time_lag_effect_mul$cr.lb, 
    ymax2 = pred_time_lag_effect_mul$cr.ub, pred = pred_time_lag_effect_mul$pred) %>% 
    ggplot(aes(x = year, y = Zr, size = (1/VZr) + 3)) + geom_point(shape = 21, 
    fill = "grey90") + # geom_ribbon(aes(ymin = ymin, ymax = ymax), fill = '#0072B2') + # not quite
# sure why this does not work
geom_smooth(aes(y = ymin2), method = "loess", se = FALSE, lty = "dotted", lwd = 0.25, 
    colour = "#0072B2") + geom_smooth(aes(y = ymax2), method = "loess", se = FALSE, 
    lty = "dotted", lwd = 0.25, colour = "#0072B2") + geom_smooth(aes(y = ymin), 
    method = "loess", se = FALSE, lty = "dotted", lwd = 0.25, colour = "#D55E00") + 
    geom_smooth(aes(y = ymax), method = "loess", se = FALSE, lty = "dotted", 
        lwd = 0.25, colour = "#D55E00") + geom_smooth(aes(y = pred), method = "loess", 
    se = FALSE, lty = "dashed", lwd = 0.5, colour = "black") + ylim(-1, 2) + 
    xlim(1994, 2019) + scale_x_continuous(breaks = c(1995, 2000, 2005, 2010, 
    2015, 2020)) + # geom_abline(intercept = mr_host_range_link_ratio$beta[[1]], slope =
# mr_host_range_link_ratio$beta[[2]], alpha = 0.7, linetype = 'dashed', size
# = 0.5) +
labs(x = "Year", y = expression(paste(italic(Zr), " (effect size)")), size = expression(paste(italic(n), 
    " (# of species)"))) + guides(fill = "none", colour = "none") + # themses
theme_bw() + theme(legend.position = c(0, 1), legend.justification = c(0, 1)) + 
    theme(legend.direction = "horizontal") + # theme(legend.background = element_rect(fill = 'white', colour = 'black'))
# +
theme(legend.background = element_blank()) + theme(axis.text.y = element_text(size = 10, 
    colour = "black", hjust = 0.5, angle = 90))

fit_time_lag_effect_mul

Supplementary Figure 10: A bubble plot showing a predicted loess line for the contentious variable year (given the values of the other 3 variables in the model) with their 95% confidences regions (orange dotted lines) and 95% prediction regions (blue dotted lines) with observed effect sizes based on various sample sizes. Note that the lines are not linear as these are based on multivariate predictions of the data points.

Sensitivity Analysis

The funnel plots above identified the issue of upper bounds of this effect size given a sample size (a upper limit of a p value given the number of randomizations). This boundary would make our estimates of mean effect sizes and contrasts (i.e., comparing two groups) so that our our overall conclusions are often conservative. To demonstrate this, we conducted two analyses to show: 1) the number of randomizations do not differ between categories in the 3 important categorical moderators (mode_of_transmission_broad, host_tax_broad, & symbiosis) and 2) categories with high effect sizes would include “bounded” effect sizes (i.e., from p = 0.01, 0.001, or 0.0001) in the 3 moderators.

Sensitivity test 1: the number of randomizations

Below, we showed that none of categorizes have significantly different numbers of randomizations in all mode_of_transmission_broad, host_tax_broad, & symbiosis.

The type of symbiosis: parasitism vs. mutualism

# 233 --- Yes = 74 (0.3175966%); No = 159

# symbiosis multiple intercepts
sa_random_symbiosis1 <- lmer(log(no_randomizations) ~ symbiosis - 1 + (1 | authors), 
    data = dat)
# contrast
sa_random_symbiosis2 <- lmer(log(no_randomizations) ~ symbiosis + (1 | authors), 
    data = dat)

Supplementary Table 19: Regression coefficients (Estimate), 95% confidence intervals (CIs), variance components (V) and variance explained, R²_[marginal] (R2) from the meta-regression with symbiont_tax_broad.

# getting marginal R2
r2_sa_random_symbiosis <- r2_nakagawa(sa_random_symbiosis1)

# getting estimates

res_sa_random_symbiosis <- tibble(estiamte = c(fixef(sa_random_symbiosis1), 
    fixef(sa_random_symbiosis2)[2]))

ci_sa_random_symbiosis1 <- confint(sa_random_symbiosis1)
ci_sa_random_symbiosis2 <- confint(sa_random_symbiosis2)
res_sa_random_symbiosis %<>% mutate(lowerCL = c(ci_sa_random_symbiosis1[3:4, 
    1], ci_sa_random_symbiosis2[4, 1]))
res_sa_random_symbiosis %<>% mutate(upperCL = c(ci_sa_random_symbiosis1[3:4, 
    2], ci_sa_random_symbiosis2[4, 2]))

# creating a table
tibble(`Fixed effect` = c(as.character(res_symbiosis1$name), cont_gen(res_symbiosis1$name)), 
    Estimate = res_sa_random_symbiosis$estiamte, `Lower CI [0.025]` = res_sa_random_symbiosis$lowerCL, 
    `Upper CI  [0.975]` = res_sa_random_symbiosis$upperCL, `V[authors]` = c(attr(VarCorr(sa_random_symbiosis1)$author, 
        "stddev")^2, rep(NA, 2)), `V[residuals]` = c(attr(VarCorr(sa_random_symbiosis1), 
        "sc")^2, rep(NA, 2)), R2 = c(r2_sa_random_symbiosis$R2_marginal, rep(NA, 
        2))) %>% kable("html", digits = 3) %>% kable_styling("striped", position = "left")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	V[authors]	V[residuals]	R2
Mutualist	7.696	7.286	8.106	2.931	0.352	0.004
Parasite	7.619	7.314	7.924	NA	NA	NA
Mutualist-Parasite	-0.077	-0.544	0.389	NA	NA	NA

The effect of host taxa

# host_tax_broad mutiple intercepts
sa_random_host_tax_broad1 <- lmer(log(no_randomizations) ~ host_tax_broad - 
    1 + (1 | authors), data = dat)
# contrast 1
sa_random_host_tax_broad2 <- lmer(log(no_randomizations) ~ host_tax_broad + 
    (1 | authors), data = dat)
# contrast 2
sa_random_host_tax_broad3 <- lmer(log(no_randomizations) ~ relevel(host_tax_broad, 
    ref = "Plant") + (1 | authors), data = dat)

# contrast 3
sa_random_host_tax_broad4 <- lmer(log(no_randomizations) ~ relevel(host_tax_broad, 
    ref = "Invert") + (1 | authors), data = dat)

Supplementary Table 20: Regression coefficients (Estimate), 95% confidence intervals (CIs), variance components (V) and variance explained, R²_[marginal] (R2) from the meta-regression with symbiont_tax_broad.

# getting marginal R2
r2_sa_random_host_tax_broad <- r2_nakagawa(sa_random_host_tax_broad1)

# getting estimates
res_sa_random_host_tax_broad <- tibble(estiamte = c(fixef(sa_random_host_tax_broad1), 
                                                    fixef(sa_random_host_tax_broad2)[2:4],
                                                    fixef(sa_random_host_tax_broad3)[3:4],
                                                    fixef(sa_random_host_tax_broad4)[4])) 
  
ci_sa_random_host_tax_broad1<-confint(sa_random_host_tax_broad1)
ci_sa_random_host_tax_broad2<-confint(sa_random_host_tax_broad2)
ci_sa_random_host_tax_broad3<-confint(sa_random_host_tax_broad3)
ci_sa_random_host_tax_broad4<-confint(sa_random_host_tax_broad4)
res_sa_random_host_tax_broad %<>% mutate(lowerCL = c(ci_sa_random_host_tax_broad1[3:6,1], 
                                                ci_sa_random_host_tax_broad2[4:6,1],
                                                ci_sa_random_host_tax_broad3[5:6,1],
                                                ci_sa_random_host_tax_broad4[6,1]))
res_sa_random_host_tax_broad %<>% mutate(upperCL = c(ci_sa_random_host_tax_broad1[3:6,2], 
                                                ci_sa_random_host_tax_broad2[4:6,2],
                                                ci_sa_random_host_tax_broad3[5:6,2],
                                                ci_sa_random_host_tax_broad4[6,2]))
# creating a table
tibble(
  `Fixed effect` =  c(as.character(res_symbiont_tax_broad1$name), cont_gen(res_symbiont_tax_broad1$name)), # done
  Estimate = res_sa_random_host_tax_broad$estiamte,
  `Lower CI [0.025]` = res_sa_random_host_tax_broad$lowerCL,
  `Upper CI  [0.975]` = res_sa_random_host_tax_broad$upperCL,
  `V[authors]` = c(attr(VarCorr(sa_random_host_tax_broad1)$author,"stddev")^2,  rep(NA, 9)),
  `V[residuals]` =c(attr(VarCorr(sa_random_host_tax_broad1),"sc")^2, rep(NA, 9)),
  `R2` = c(r2_sa_random_host_tax_broad$R2_marginal, rep(NA, 9))) %>% kable("html", digits = 3) %>%
  kable_styling("striped", position = "left") 

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	V[authors]	V[residuals]	R2
Microbe	8.223	7.279	9.168	2.937	0.35	0.071
Plant	7.466	6.856	8.076	NA	NA	NA
Invert	7.650	7.167	8.134	NA	NA	NA
Vert	7.634	7.265	8.003	NA	NA	NA
Microbe-Plant	-0.757	-1.882	0.367	NA	NA	NA
Microbe-Invert	-0.573	-1.634	0.488	NA	NA	NA
Microbe-Vert	-0.589	-1.603	0.425	NA	NA	NA
Plant-Invert	0.184	-0.594	0.962	NA	NA	NA
Plant-Vert	0.168	-0.545	0.881	NA	NA	NA
Invert-Vert	-0.016	-0.599	0.566	NA	NA	NA

The effect of the model of transmission

# mode_of_transmission_broad

sa_random_mode_of_transmission_broad1 <- lmer(log(no_randomizations) ~ mode_of_transmission_broad - 
    1 + (1 | authors), data = dat)


# contrast 1
sa_random_mode_of_transmission_broad2 <- lmer(log(no_randomizations) ~ mode_of_transmission_broad + 
    (1 | authors), data = dat)

# contrast 2
sa_random_mode_of_transmission_broad3 <- lmer(log(no_randomizations) ~ relevel(mode_of_transmission_broad, 
    ref = "vertical") + (1 | authors), data = dat)

Supplementary Table 21: Regression coefficients (Estimate), 95% confidence intervals (CIs), variance components (V) and variance explained, R²_[marginal] (R2) from the meta-regression with symbiont_tax_broad.

# getting marginal R2
r2_sa_random_mode_of_transmission_broad <- r2_nakagawa(sa_random_mode_of_transmission_broad1)

# getting estimates
res_sa_random_mode_of_transmission_broad <- tibble(estiamte = c(fixef(sa_random_mode_of_transmission_broad1), 
    fixef(sa_random_mode_of_transmission_broad2)[2:3], fixef(sa_random_mode_of_transmission_broad3)[3]))

ci_sa_random_mode_of_transmission_broad1 <- confint(sa_random_mode_of_transmission_broad1)
ci_sa_random_mode_of_transmission_broad2 <- confint(sa_random_mode_of_transmission_broad2)
ci_sa_random_mode_of_transmission_broad3 <- confint(sa_random_mode_of_transmission_broad3)
res_sa_random_mode_of_transmission_broad %<>% mutate(lowerCL = c(ci_sa_random_mode_of_transmission_broad1[3:5, 
    1], ci_sa_random_mode_of_transmission_broad2[4:5, 1], ci_sa_random_mode_of_transmission_broad3[5, 
    1]))
res_sa_random_mode_of_transmission_broad %<>% mutate(upperCL = c(ci_sa_random_mode_of_transmission_broad1[3:5, 
    2], ci_sa_random_mode_of_transmission_broad2[4:5, 2], ci_sa_random_mode_of_transmission_broad3[5, 
    2]))
# creating a table
tibble(`Fixed effect` = c(as.character(res_mode_of_transmission_broad1$name), 
    cont_gen(res_mode_of_transmission_broad1$name)), Estimate = res_sa_random_mode_of_transmission_broad$estiamte, 
    `Lower CI [0.025]` = res_sa_random_mode_of_transmission_broad$lowerCL, `Upper CI  [0.975]` = res_sa_random_mode_of_transmission_broad$upperCL, 
    `V[authors]` = c(attr(VarCorr(sa_random_mode_of_transmission_broad1)$author, 
        "stddev")^2, rep(NA, 5)), `V[residuals]` = c(attr(VarCorr(sa_random_mode_of_transmission_broad1), 
        "sc")^2, rep(NA, 5)), R2 = c(r2_sa_random_mode_of_transmission_broad$R2_marginal, 
        rep(NA, 5))) %>% kable("html", digits = 3) %>% kable_styling("striped", 
    position = "left")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	V[authors]	V[residuals]	R2
both	7.535	6.978	8.091	2.913	0.363	0.12
horizontal	7.515	7.164	7.866	NA	NA	NA
vertical	8.089	7.509	8.668	NA	NA	NA
both-horizontal	-0.020	-0.678	0.638	NA	NA	NA
both-vertical	0.554	-0.250	1.357	NA	NA	NA
horizontal-vertical	-0.574	-1.251	0.104	NA	NA	NA

Sensitivity test 2: reaching the limits

Below, we showed that categorizes with higher effect sizes were more likely to have “bounded” effect sizes (limit_rearched) in all mode_of_transmission_broad, host_tax_broad, & symbiosis. This indicate that the higher the estimate of mean effect size is, the more underestimated the mean effect size is. This is true for differences between two categories; the larger the difference between two, the more underestimated the difference is.

The type of symbiosis: parasitism vs. mutualism

Supplementary Table 22: Regression coefficients (Estimate), 95% confidence intervals (CIs), variance components (V) and variance explained, R²_[marginal] (R2) from the meta-regression with symbiont_tax_broad.

# symbiosis
sa_limit_symbiosis1 <- glmer(limit_rearched ~ symbiosis - 1 + (1 | authors), 
    family = "binomial", data = dat)

# getting marginal R2
r2_sa_limit_symbiosis <- r2_nakagawa(sa_limit_symbiosis1)

# getting estimates
res_sa_limit_symbiosis <- tibble(estiamte = fixef(sa_limit_symbiosis1))

res_sa_limit_symbiosis %<>% mutate(lowerCL = (tidy(sa_limit_symbiosis1)$estimate[-3] - 
    tidy(sa_limit_symbiosis1)$std.error[-3] * qnorm(0.975)))
res_sa_limit_symbiosis %<>% mutate(upperCL = (tidy(sa_limit_symbiosis1)$estimate[-3] + 
    tidy(sa_limit_symbiosis1)$std.error[-3] * qnorm(0.975)))
# creating a table
tibble(`Fixed effect` = as.character(res_symbiosis1$name), Estimate = res_sa_limit_symbiosis$estiamte, 
    `Lower CI [0.025]` = res_sa_limit_symbiosis$lowerCL, `Upper CI  [0.975]` = res_sa_limit_symbiosis$upperCL, 
    `V[authors]` = c(attr(VarCorr(sa_limit_symbiosis1)$author, "stddev")^2, 
        rep(NA, 1)), R2 = c(r2_sa_limit_symbiosis$R2_marginal, rep(NA, 1))) %>% 
    kable("html", digits = 3) %>% kable_styling("striped", position = "left")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	V[authors]	R2
Mutualist	-0.401	-1.074	0.272	1.657	0.051
Parasite	-1.309	-2.015	-0.603	NA	NA

The effect of host taxa

Supplementary Table 23: Regression coefficients (Estimate), 95% confidence intervals (CIs), variance components (V) and variance explained, R²_[marginal] (R2) from the meta-regression with symbiont_tax_broad.

# host_tax_broad
sa_limit_host_tax_broad1 <- glmer(limit_rearched ~ host_tax_broad - 1 + (1 | 
    authors), family = "binomial", data = dat)

# getting marginal R2
r2_sa_limit_host_tax_broad <- r2_nakagawa(sa_limit_host_tax_broad1)

# getting estimates
res_sa_limit_host_tax_broad <- tibble(estiamte = fixef(sa_limit_host_tax_broad1))

res_sa_limit_host_tax_broad %<>% mutate(lowerCL = (tidy(sa_limit_host_tax_broad1)$estimate[-5] - 
    tidy(sa_limit_host_tax_broad1)$std.error[-5] * qnorm(0.975)))
res_sa_limit_host_tax_broad %<>% mutate(upperCL = (tidy(sa_limit_host_tax_broad1)$estimate[-5] + 
    tidy(sa_limit_host_tax_broad1)$std.error[-5] * qnorm(0.975)))
# creating a table
tibble(`Fixed effect` = as.character(res_symbiont_tax_broad1$name), Estimate = res_sa_limit_host_tax_broad$estiamte, 
    `Lower CI [0.025]` = res_sa_limit_host_tax_broad$lowerCL, `Upper CI  [0.975]` = res_sa_limit_host_tax_broad$upperCL, 
    `V[authors]` = c(attr(VarCorr(sa_limit_host_tax_broad1)$author, "stddev")^2, 
        rep(NA, 3)), R2 = c(r2_sa_limit_host_tax_broad$R2_marginal, rep(NA, 
        3))) %>% kable("html", digits = 3) %>% kable_styling("striped", position = "left")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	V[authors]	R2
Microbe	-1.211	-2.737	0.316	1.358	0.035
Plant	-1.574	-2.571	-0.578	NA	NA
Invert	-0.579	-1.308	0.150	NA	NA
Vert	-0.937	-1.588	-0.286	NA	NA

The effect of the model of transmission

Supplementary Table 24: Regression coefficients (Estimate), 95% confidence intervals (CIs), variance components (V) and variance explained, R²_[marginal] (R2) from the meta-regression with symbiont_tax_broad.

# mode_of_transmission_broad
sa_limit_mode_of_transmission_broad1 <- glmer(limit_rearched ~ mode_of_transmission_broad - 
    1 + (1 | authors), family = "binomial", data = dat)

# getting marginal R2
r2_sa_limit_mode_of_transmission_broad <- r2_nakagawa(sa_limit_mode_of_transmission_broad1)

# getting estimates
res_sa_limit_mode_of_transmission_broad <- tibble(estiamte = fixef(sa_limit_mode_of_transmission_broad1))

res_sa_limit_mode_of_transmission_broad %<>% mutate(lowerCL = (tidy(sa_limit_mode_of_transmission_broad1)$estimate[-4] - 
    tidy(sa_limit_mode_of_transmission_broad1)$std.error[-4] * qnorm(0.975)))
res_sa_limit_mode_of_transmission_broad %<>% mutate(upperCL = (tidy(sa_limit_mode_of_transmission_broad1)$estimate[-4] + 
    tidy(sa_limit_mode_of_transmission_broad1)$std.error[-4] * qnorm(0.975)))
# creating a table
tibble(`Fixed effect` = as.character(res_mode_of_transmission_broad1$name), 
    Estimate = res_sa_limit_mode_of_transmission_broad$estiamte, `Lower CI [0.025]` = res_sa_limit_mode_of_transmission_broad$lowerCL, 
    `Upper CI  [0.975]` = res_sa_limit_mode_of_transmission_broad$upperCL, `V[authors]` = c(attr(VarCorr(sa_limit_mode_of_transmission_broad1)$author, 
        "stddev")^2, rep(NA, 2)), R2 = c(r2_sa_limit_mode_of_transmission_broad$R2_marginal, 
        rep(NA, 2))) %>% kable("html", digits = 3) %>% kable_styling("striped", 
    position = "left")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	V[authors]	R2
both	-0.928	-1.809	-0.046	1.427	0.078
horizontal	-1.326	-2.031	-0.621	NA	NA
vertical	0.060	-0.750	0.869	NA	NA

Recommendations based on our sensitivity analysis

We have 3 recommendations for future (and past) co-divergence work using TreeMap and ParaFit.

• Researchers need to set the number of randomizations much higher to obtain more accurate p value.
• If possible, the reanalysis of earlier trees are necessary; also, analysis using updated trees will be desirable.
• New studies should provide all data available online so that a future meta-analysis can incorporate all available data.

Acknowledgement

Many coding materials have been borrowed from these paper^14,15. We thank Losia Lagisz for preparing small icons and cartoons used in the figures.

R Session Information

sessionInfo() %>% pander()

R version 3.6.1 (2019-07-05)

Platform: x86_64-apple-darwin15.6.0 (64-bit)

locale: en_AU.UTF-8||en_AU.UTF-8||en_AU.UTF-8||C||en_AU.UTF-8||en_AU.UTF-8

attached base packages: grid, stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: patchwork(v.0.0.1), png(v.0.1-7), performance(v.0.3.0), broom.mixed(v.0.2.4), lme4(v.1.1-21), MuMIn(v.1.43.6), plotly(v.4.9.0), ggbeeswarm(v.0.6.0), MCMCglmm(v.2.29), ape(v.5.3), coda(v.0.19-3), metafor(v.2.1-0), Matrix(v.1.2-17), pander(v.0.6.3), magrittr(v.1.5), gridExtra(v.2.3), kableExtra(v.1.1.0), forcats(v.0.4.0), stringr(v.1.4.0), dplyr(v.0.8.3), purrr(v.0.3.2), readr(v.1.3.1), tidyr(v.0.8.3), tibble(v.2.1.3), ggplot2(v.3.2.1) and tidyverse(v.1.2.1)

loaded via a namespace (and not attached): httr(v.1.4.1), jsonlite(v.1.6), viridisLite(v.0.3.0), splines(v.3.6.1), modelr(v.0.1.5), assertthat(v.0.2.1), highr(v.0.8), stats4(v.3.6.1), tensorA(v.0.36.1), vipor(v.0.4.5), cellranger(v.1.1.0), bayestestR(v.0.2.5), yaml(v.2.2.0), pillar(v.1.4.2), backports(v.1.1.4), lattice(v.0.20-38), glue(v.1.3.1), digest(v.0.6.20), rvest(v.0.3.4), minqa(v.1.2.4), colorspace(v.1.4-1), plyr(v.1.8.4), htmltools(v.0.3.6), pkgconfig(v.2.0.2), broom(v.0.5.2), haven(v.2.1.1), corpcor(v.1.6.9), scales(v.1.0.0), webshot(v.0.5.1), cubature(v.2.0.3), generics(v.0.0.2), pacman(v.0.5.1), withr(v.2.1.2), TMB(v.1.7.15), lazyeval(v.0.2.2), cli(v.1.1.0), crayon(v.1.3.4), readxl(v.1.3.1), evaluate(v.0.14), nlme(v.3.1-141), MASS(v.7.3-51.4), xml2(v.1.2.2), beeswarm(v.0.2.3), tools(v.3.6.1), data.table(v.1.12.2), hms(v.0.5.0), formatR(v.1.7), munsell(v.0.5.0), compiler(v.3.6.1), rlang(v.0.4.0), nloptr(v.1.2.1), rstudioapi(v.0.10.0-9000), htmlwidgets(v.1.3), labeling(v.0.3), base64enc(v.0.1-3), rmarkdown(v.1.14), boot(v.1.3-23), codetools(v.0.2-16), gtable(v.0.3.0), reshape2(v.1.4.3), R6(v.2.4.0), lubridate(v.1.7.4), knitr(v.1.24), zeallot(v.0.1.0), insight(v.0.4.1), stringi(v.1.4.3), parallel(v.3.6.1), Rcpp(v.1.0.2), vctrs(v.0.2.0), tidyselect(v.0.2.5) and xfun(v.0.8)

References

1. Page, R. D. Parallel phylogenies: Reconstructing the history of host-parasite assemblages. Cladistics 10, 155–173 (1994).

2. Legendre, P., Desdevises, Y. & Bazin, E. A statistical test for host–parasite coevolution. Systematic biology 51, 217–234 (2002).

3. Rosenthal, R. & Rubin, D. B. Requivalent: A simple effect size indicator. Psychological methods 8, 492 (2003).

4. Viechtbauer, W. Conducting meta-analyses in R with the metafor package. Journal of Statistical Software 36, 1–48 (2010).

5. Knapp, G. & Hartung, J. Improved tests for a random effects meta-regression with a single covariate. Statistics in medicine 22, 2693–2710 (2003).

6. Higgins, J. P., Thompson, S. G., Deeks, J. J. & Altman, D. G. Measuring inconsistency in meta-analyses. Bmj 327, 557–560 (2003).

7. Nakagawa, S. & Santos, E. S. Methodological issues and advances in biological meta-analysis. Evolutionary Ecology 26, 1253–1274 (2012).

8. Senior, A. M. et al. Heterogeneity in ecological and evolutionary meta-analyses: Its magnitude and implications. Ecology 97, 3293–3299 (2016).

9. Nakagawa, S. & Schielzeth, H. A general and simple method for obtaining r2 from generalized linear mixed-effects models. Methods in ecology and evolution 4, 133–142 (2013).

10. Barton, K. MuMIn: Multi-model inference. http://r-forge. r-project. org/projects/mumin/ (2009).

11. Peters, J. L., Sutton, A. J., Jones, D. R., Abrams, K. R. & Rushton, L. Contour-enhanced meta-analysis funnel plots help distinguish publication bias from other causes of asymmetry. Journal of clinical epidemiology 61, 991–996 (2008).

12. Egger, M., Smith, G. D., Schneider, M. & Minder, C. Bias in meta-analysis detected by a simple, graphical test. Bmj 315, 629–634 (1997).

13. Moreno, S. G. et al. Assessment of regression-based methods to adjust for publication bias through a comprehensive simulation study. BMC medical research methodology 9, 2 (2009).

14. Cally, J. G., Stuart-Fox, D. & Holman, L. Meta-analytic evidence that sexual selection improves population fitness. Nature communications 10, 2017 (2019).

15. O’Dea, R. E., Lagisz, M., Hendry, A. P. & Nakagawa, S. Developmental temperature affects phenotypic means and variability: A meta-analysis of fish data. Fish and Fisheries doi.org/10.1111/faf.12394, (2019).